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Abstract. Investigations of emergent symmetry breaking phenomena occurring in 

small finite-size systems are reviewed, with a focus on the strongly correlated regime 
of electrons in two-dimensional semicoductor quantum dots and trapped ultracold 
bosonic atoms in harmonic traps. Throughout the review we emphasize universal 
aspects and similarities of syninietry breaking found in these systems, as well as in more 
traditional fields like nuclear physics and quantum chemistry, which are characterized 
by very different interparticle forces. A unified description of strongly correlated 
phenomena in finite systems of repelling particles (whether fermions or bosons) is 
presented through the development of a two-step method of symmetry breaking at 
the unrestricted Hartree-Fock level and of subsequent symmetry restoration via post 
Hartree-Fock projection techniques. Quantitative and qualitative aspects of the two- 
step method are treated and validated by exact diagonalization calculations. 
Strongly-correlated phenomena emerging from symmetry breaking include: 

(I) Chemical bonding, dissociation, and entanglement (at zero and finite magnetic 
fields) in quantum dot molecules and in pinned electron molecular dimers formed 
within a single anisotropic quantum dot, with potential technological applications to 
solid-state quantum-computing devices. 

(II) Electron crystallization, with particle localization on the vertices of concentric 
polygonal rings, and formation of rotating electron molecules (REMs) in circular 
quantum dots. Such electron molecules exhibit ro- vibrational excitation spectra, in 
analogy with natural molecules. 

(III) At high magnetic fields, the REMs are described by parameter-free analytic wave 
functions, which are an alternative to the Laughlin and composite-fermion approaches, 
offering a new point of view of the fractional quantum Hall regime in quantum dots 
(with possible implications for the thermodynamic limit). 

(IV) Crystalline phases of strongly repelling bosons. In rotating traps and in analogy 

with the REMs, such repelling bosons form rotating boson molecules (RBMs). For a 
small number of bosons, the RBMs are energetically favored compared to the Gross- 
Pitaevskii solutions describing vortex formation. 

We discuss the present status concerning experimental signatures of such strongly 
correlated states, in view of the promising outlook created by the latest experimental 
improvements that are achieving unprecedented control over the range and strength of 
interparticle interactions. 
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1. Introduction 

1.1. Preamble 

Fermionic or bosonic particles confined in manmade devices, i.e., electrons in two- 
dimensional (2D) quantum dots (QDs, referred to also as artificial atoms) or ultracold 
atoms in harmonic traps, can localize and form structures with molecular, or crystalline, 
characteristics. These molecular states of localized particles differ in an essential 
way from the electronic-shell-structure picture of delocalized electrons filling successive 
orbitals in a central-mean-field potential (the Aufbau principle), familiar from the many- 
body theory of natural atoms and the Mendeleev periodic table; they also present a 
different regime from that exhibited by a Bose-Einstein condensate (BEC, associated 
often with the mean-field Gross-Pitaevskii equation). The molecular states originate 
from strong correlations between the constituent repelling particles and they are called 
electron (and often Wigner) or boson molecules. 

Such molecular states forming within a single confining potential well constitute new 
phases of matter and allow for investigations of novel strongly-correlated phenomena 
arising in physical systems with a range of materials' characteristics unavailable 
experimentally (and theoretically unexplored) until recently. One example is the range 
of values of the socalled Wigner parameter (denoted as R\y for charged particles and Rs 



for neutral ones, see section 2.1.2) which expresses the relative strength of the two-body 
repulsion and the one-particle kinetic energy, reflecting and providing a measure of the 
strength of correlations in the system under study. For the two-dimensional systems 
which we discuss here, these values are often larger than the corresponding ones for 
natural atoms and molecules. 

Other research opportunities offered by the quantum-dot systems are related to 
their relatively large (spatial) size (arising from a small electron effective mass and 
large dielectric constant), which allows the full range of orbital magnetic effects to be 
covered for magnetic fields that are readily attained in the laboratory (less than 40 T). 
In contrast, for natural atoms and molecules, magnetic fields of sufficient strength (i.e., 
larger than 10^ T) to produce novel phenomena related to orbital magnetism (beyond 
the perturbative regime) are known to occur only in astrophysical environments (e.g., 
on the surface of neutron stars) [1] . For ultracold gases, a similar extraordinary physical 
regime can be reached via the fast rotation of the harmonic trap. 

In addition to the fundamental issues unveiled through investigations of molecular 
states in quantum dots, these strongly-correlated states are of technological significance 
because of the potential use of manmade nanoscale systems for the implementation of 
qubits and quantum logic gates in quantum computers. 

The existence of electron and boson molecules is supported by large-scale exact 
diagonalization (EXD) calculations, which provide the ultimate theoretical test. The 
discovery of these "crystalline" states has raised important fundamental aspects, 
including the nature of quantum phase transitions and the conceptual issues relating to 
spontaneous symmetry breaking (SSB) in small finite-size systems. 
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The present report addresses primarily the physics, theoretical description, and 
fundamental many-body aspects of molecular (crystalline) states in small systems. For 
a comprehensive description of the electronic-shell-structure regime (Aufbau-principle 
regime) in quantum dots and of Bose-Einstein condensates in harmonic traps, see the 
earlier reviews by Kouwenhoven et al [2] (QDs), Reimann and Manninen [3] (QDs), 
Dalfovo et al [1] (BECs), and Leggett [5] (BECs). Furthermore, in larger quantum 
dots, the symmetries of the external confinement that lead to shell structure are broken, 
and such dots exhibit mesoscopic fluctuations and interplay between single-particle 
quantum chaos [6] and many-body correlations. For a comprehensive description of 
this mesoscopic regime in quantum dots, see the reviews by Beenakker [7] and Alhassid 
0. 

1.2. Spontaneous symmetry breaking: confined geometries versus extended systems 

Spontaneous symmetry breaking is a ubiquitous phenomenon in the macroscopic 
world. Indeed, there is an abundance of macroscopic systems and objects that are 
observed, or can be experimentally prepared, with effective many-body ground states 
whose symmetry is lower than the symmetry of the underlying many-body quantum- 
mechanical Hamiltonian; one says that in such cases the system lowers its energy through 
spontaneous symmetry breaking, resulting in a state of lower symmetry and higher order. 
It is important to stress that macroscopic SSB strongly suppresses quantum fluctuations 
and thus it can be described appropriately by a set of non-linear mean-field equations 
for the "order parameter." The appearance of the order parameter is governed by 
bifurcations associated with the non-linearity of the mean-field equations and has led to 
the notion of "emergent phenomena," a notion that helped promote condensed-matter 
physics as a branch of physics on a par with high-energy particle physics (in reference to 
the fundamental nature of the pursuit in these fields; see the seminal paper by Anderson 
in Ref. |9]). 

Our current understanding of the physics of SSB in the thermodynamic limit (when 
the number of particles N — > oo) owes a great deal to the work of Anderson [10], who 
suggested that the broken-symmetry state can be safely taken as the effective ground 
state. In arriving at this conclusion Anderson invoked the concept of (generalized) 
rigidity. As a concrete example, one would expect a crystal to behave like a macroscopic 
body, whose Hamiltonian is that of a heavy rigid rotor with a low-energy excitation 
spectrum L'^/2J' of angular-momentum (L) eigenstates, with the moment of inertia 
being of order (macroscopically large when N — > oo). The low-energy excitation 
spectrum of this heavy rigid rotor above the ground-state (L = 0) is essentially gapless 
(i.e., continuous). Thus although the formal ground state posseses continuous rotational 
symmetry (i.e., L = 0), "there is a manifold of other states, degenerate in the N ^ oo 
limit, which can be recombined to give a very stable wave packet with essentially the 
nature" of the broken-symmetry state (see p 44 in Ref. jTU]). 

As a consequence of the "macroscopic heaviness" as ^ oo, the relaxation of 
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the system from the wave packet state (i.e., the broken-symmetry state) to the exact 
symmetrical ground state becomes exceedingly long. Consequently, in this limit, when 
symmetry breaking occurs, there is practically no need to follow up with a symmetry 
restoration step; that is the symmetry-broken state is admissible as an effective ground 
state. 

The present report addresses the much less explored question of symmetry breaking 
in finite condensed-matter systems with a small number of particles. For small systems, 
spontaneous symmetry breaking appears again at the level of mean-field description 
[e.g., the Hartree-Fock (HF) level]. A major difference from the N ^ oo limit, however, 
arises from the fact that quantum fluctuations in small systems cannot be neglected. 
To account for the large fluctuations, one has to perform a subsequent post-Hartree- 
Fock step that restores the broken symmetries (and the linearity of the many-body 
Schrodinger equation). Subsequent to symmetry restoration, the ground state obeys all 
the original symmetries of the many-body Hamiltonian; however, effects of the mean- 
field symmetry breaking do survive in the properties of the ground state of small systems 
and lead to emergent phenomena associated with formation of novel states of matter 
and with characteristic behavior in the excitation spectra. In the following, we will 
present an overview of the current understanding of SSB in small systems focusing 
on the essential theoretical aspects, as well as on the contributions made by SSB-based 
approaches to the fast developing fields of two-dimensional semiconductor quantum dots 
and ultracold atomic gases in harmonic and toroidal traps. 

1.3. Historical background from nuclear physics and chemistry 

The mean field approach, in the form of the Hartree-Fock theory and of the Gross- 
Pitaevskii (GP) equation, has been a useful tool in elucidating the physics of finite- 
size fermionic and bosonic systems, respectively. Its applications cover a wide range 
of systems, from natural atoms, natural molecules, and atomic nuclei, to metallic 
nanoclusters, and most recently two-dimensional quantum dots and ultracold gases 
confined in harmonic (parabolic) traps. Of particular interest for the present review (due 
to spatial-symmetry-breaking aspects) has been the mean-field description of deformed 
nuclei [HI |12l [13] and metal clusters [TH HH HE] (exhibiting ellipsoidal shapes). At a 
first level of description, deformation effects in these latter systems can be investigated 
via semi-empirical mean-field models, like the particle-rotor model [H] of Bohr and 
Mottelson (nuclei), the anisotropic-harmonic-oscillator model of Nilsson (nuclei [12] and 
metal clusters [S]), and the shell-correction method of Strutinsky (nuclei [T7] and metal 
clusters [T5| IT6]). At the microscopic level, the mean field for fermions is often described 
[IS] [TU] via the self-consistent single-determinantal Hartree-Fock theory. At this level, 
the description of deformation effects mentioned above requires [H] consideration of 
unrestricted Hartree-Fock (UHF) wave functions that break explicitly the rotational 
symmetries of the original many-body Hamiltonian, but yield HF Slater determinants 
with lower energy compared to the symmetry-adapted restricted Hartree-Fock (RHF) 
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solutions It] 

In earlier publications [201 1211 1221 l23l l2ll ISSl |26] , we have shown that, in the strongly 
correlated regime, UHF solutions that violate the rotational (circular) symmetry arise 
most naturally in the case of two-dimensional single quantum dots, for both the cases 
of zero and high magnetic field; for a UHF calculation in the lowest Landau level 
(LLL), see also Ref. [21] • Unlike the case of atomic nuclei, however, where (due to 
the attractive interaction) symmetry breaking is associated primarily with quadrupole 
shape deformations (a type of Jahn- Teller distortion), spontaneous symmetry breaking 
in 2D quantum dots induces electron localization (or "crystallization") associated with 
formation of electron, or Wigner, molecules). The latter name is used in honor of 
Eugene Wigner who predicted the formation of a classical rigid Wigner crystal for the 
3D electron gas at very low densities [28j. We stress, however, that because of the 
finite size, Wigner molecules are most often expected to show a physical behavior quite 
different from the classical Wigner crystal. Indeed, for finite N, Wigner molecules 
exhibit analogies closer to natural molecules, and the Wigner-crystal limit is expected 
to be reached only for special limiting conditions. 

For a small system the violation in the mean-field approximation of the symmetries 
of the original many-body Hamiltonian appears to be paradoxical at a first glance, 
and some times it has been described mistakenly as an "artifact" (in particular in the 
context of density- functional theory [29]). However, for the specific cases arising in 
Nuclear Physics and Quantum Chemistry, two theoretical developments had already 
resolved this paradox. They are: (1) the theory of restoration of broken symmetries via 
projection technique^ [301 ED 132], and (2) the group theoretical analysis of symmetry- 
broken HF orbitals and solutions in chemical reactions, initiated by Fukutome and 
coworkers [33] who used the symmetry groups associated with the natural 3D molecules. 
Despite the different fields, the general principles established in these earlier theoretical 
developments in nuclear physics and quantum chemistry have provided a wellspring of 
assistance in our investigations of symmetry breaking for electrons in quantum dots and 
bosons in harmonic traps. In particular, the restoration of broken symmetries in QDs 
and ultracold atomic traps via projection techniques constitutes a main theme of the 
present report. 

The theory of restoration of broken symmetries has been developed into a 
sophisticated computational approach in modern nuclear physics. Using the broken- 
symmetry solutions of the Hartree-Fock-Bogoliubov theorj|§] (that accounts for nuclear 
pairing and superfluidity), this approach has been proven particularly efficient in 
describing the competition between shape deformation and pairing in nuclei. For some 
recent papers in nuclear physics, see, e.g., Refs. [3H |35l |36l [371 [381 139] ; for application 

t See in particular Ch 5.5 and Ch 11 in Ref. [TS]. However, our terminology (i.e., UHF vs. RHF) 
follows the practice in quantum chemistry (see Ref. [H]). 

I For the restoration of broken rotational symmetries in atomic nuclei, see Ref. |3D| and Ch 1 1 in Ref. 
[TH] . For the restoration of broken spin symmetries in natural 3D molecules, see Ref. [3T] . 
§ See Ch 7 in Ref. [TS;. 
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to superconducting metallic grains, see Ref. [ID]. Pairing effects arise only in the case 
of attractive interactions and they are not considered in this report, since we deal only 
with repulsive two-body interactions. 

1.4- Scope of the review 

Having discussed earlier the general context and historical background from other fields 
regarding symmetry breaking, we give here an outline of the related methodologies and 
of the newly discovered strongly correlated phenomena that are discussed in this report 
in the area of condensed-matter nanosystems. 

In particular, a two-step method [20l [211 1221 [231 [Ml [25] of symmetry breaking at the 
unrestricted Hartree-Fock level and of subsequent post-Hartree-Fock restoration of the 
broken symmetries via projection techniques is reviewed for the case of two-dimensional 
(2D) semiconductor quantum dots and ultracold bosons in rotating traps with a small 
number (N) of particles. The general principles of the two-step method can be traced 
to nuclear theory (Peierls and Yoccoz, see the original Ref. [30], but also the recent Refs. 
[311[3S1[M1[SZ1[SSII3S]) and quantum chemistry (Lowdin, see Ref. [SI]); in the context of 
condensed-matter nanophysics and the physics of ultracold atomic gases, it constitutes a 
novel powerful many-body approach that has led to unexpected discoveries in the area of 
strongly correlated phenomena. The successes of the method have generated a promising 
theoretical outlook, bolstered by the unprecedented experimental and technological 
advances, pertaining particularly to control of system parameters (most importantly 
of the strength and variety of two-body interactions), that can be achieved in manmade 
nanostructures. 

In conjunction with exact diagonalization calculations [261 [ID [121 [131 [S] and recent 
experiments [HI [HI [33], it is shown that the two-step method can describe a wealth 
of novel strongly correlated phenomena in quantum dots and ultracold atomic traps. 
These include: 

(I) Chemical bonding, dissociation, and entanglement in quantum dot molecules 
[20] [22] [l6] and in electron molecular dimers formed within a single elliptic QD 
[m [121 [131 HI] , with potential technological applications to solid-state quantum logic 
gates [3TI[1H1I19]. 

(II) Electron crystallization, with localization on the vertices of concentric polygonal 
rings, and formation of rotating electron molecules (REMs) in circular QDs. At zero 
magnetic field (B), the REMs can approach the limit of a rigid rotor [50l |5T] : at high 
B, the REMs are highly floppy and "supersolid"-like, that is, they exhibit [SH [321 ES] 
a non-rigid rotational inertia [51], with the rings rotating independently of each other 
[321 [33]. 

(III) At high magnetic fields and under the restriction of the many-body Hilbert 
space to the lowest Landau level, the two-step method yields fully analytic many-body 
wave functions [211 [2S], which are an alternative to the Jastrow/Laughlin (JL) [33] and 
composite- fermion (CF) [561 [37] approaches, offering a new point of view of the fractional 
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quantum Hall regime (FQHE) [5H1 EH] in quantum dots (with possible implications for 
the thermodynamic limit). 

Large scale exact-diagonalization calculations [261 ESI EH] support the results of the 
two-step method outlined in items II and III above. 

(IV) The two-step method has been used jHO] to discover crystalline phases of 
strongly repelling ultracold bosons (impenetrable bosons/ Tonks- Girardeau regime 
[UT| in2]) in 2D harmonic traps. In the case of rotating traps, such repelling bosons 
form rotating boson molecules (RBMs) |63] that are energetically favorable compared 
to the Gross-Pitaevkii solutions, even for weak repulsion and, in particular, in the regime 
of GP vortex formation. 

We will not discuss in this report specific applications of the two-step method to 
atomic nuclei. Rather, as the title conveys, the report aims at exploring the universal 
characteristics of quantum correlations arising from symmetry breaking across various 
fields dealing with small finite systems, such as 2D quantum dots, trapped ultracold 
atoms, and nuclei - and even natural 3D molecules. Such universal characteristics and 
similarities in related methodologies persist across the aforementioned fields in spite of 
the differences in the size of the physical systems and in the range, nature, and strength 
of the two-body interactions. For specific applications to atomic nuclei, the interested 
reader is invited to consult the nucler physics literature cited in this report. 

1.5. Using a hierarchy of approximations versus probing of exact solutions 

Figure [T] presents a synopsis of the hierarchy of approximations associated with the two- 
step method, and in particular for the case of 2D quantum dots. (A similar synopsis can 
also be written for the case of bosonic systems.) This method produces approximate 
wave functions with lower energy at each approximation level (as indicated by the 
downward vertical arrow on the left of the figure). 

At the lowest level of approximation (corresponding to higher energy with no 
correlations included), one places the restricted Hartree-Fock, whose main restriction 
is the double occupancy (up and down spins) of each space orbital. The many-body 
wave function is a single Slater determinant associated with a "central mean field." The 
RHF preserves all spin and space symmetries. For 2D quantum dots, the single-particle 
density [also referred to as electron density (e-density)] is circularly symmetric. 

The next approximation involves the unrestricted Hartree-Fock, which employs 
different space orbitals for the two different spin directions. The UHF preserves the spin 
projection, but allows the total-spin and space symmetries (i.e., rotational symmetries 
or parity) to be broken. The broken symmetry solutions, however, are not devoid of any 
symmetry; they exhibit characteristic lower symmetries (point-group symmetries) that 
are explicit in the electron densities. The UHF many-body wave function is a single 
Slater determinant associated with a "non-central mean field." 

Subsequent approximations aim at restoring the broken symmetries via projection 
techniques. The restoration-of-symmetry step goes beyond the mean field approximation 
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Restoration of linearity 
of many-body equatons 



Figure 1. (Color online) Synopsis of the method of hierarchical approximations 
(also referred to as the "two-step method," emphasizing that symmetry breaking at 
the mean-field level must be accompanied by a subsequent post-Hartree-Fock step of 
symmetry restoration, with a subsequent further lowering of the energy). See text for 
a detailed description. 



and it provides a many-body wave function |$^^'^) that is a linear superposition of 
Slater determinants (see detailed description in section 2.2 below). The projected (PRJ) 
many-body wave function 1$^^'') preserves all the symmetries of the original many-body 
Hamiltonian; it has good total spin and angular momentum quantum numbers, and as 
a result the circular symmetry of the electron densities is restored. 

However, the lower (point-group) spatial symmetry found at the broken-symmetry 
UHF level (corresponding to the first step in this method) does not disappear. Instead, 
it becomes intrinsic or hidden, and it can be revealed via an inspection of conditional 
probability distributions (CPDs), defined as (within a proportionality constant) 

P(r, ro) = (<I>P^J| E ^(r. " r)5(r, - ro)!*^^^), (1.1) 

where $™'^(ri, r2, . . . , r^v) denotes the projected many-body wave function under 
consideration. 

If one needs to probe the intrinsic spin distribution of the localized electrons, one has 
to consider spin-resolved two-point correlation functions (spin- resolved CPDs), defined 
as 

P..„(r,ro) = (<l>™J|E5(r-r,)5(ro-r,)5..,^5.„.J$Pi^J). (1.2) 
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The spin-resolved CPD gives the spatial probability distribution of finding a second 
electron with spin projection cr under the condition that a first electron is located (fixed) 
at ro with spin projection do; a and ctq can be either up (j) or down (|). The meaning 
of the space-only CPD in (1.1) is analogous, but without consideration of spin. 

Further signatures of the intrinsic lower symmetry occur in the excitation spectra 
of circular quantum dots that exhibit ro-vibrational character related to the intrinsic 
molecular structure, or in the dissociation of quantum dot molecules. 

As the scheme in figure [T] indicates, the mean- field HF equations are non-linear and 
the symmetry breaking is associated with the appearance of bifurcations in the total HF 
energies. The occurrence of such bifurcations cannot be predicted a priori from a mere 
inspection of the many-body Hamiltonian itself; it is a genuine many-body effect that 
belongs to the class of so-called emergent phenomena P, EH ES] that may be revealed 
only through the solutions of the Hamiltonian themselves (if obtainable) or through 
experimental signatures. We note that the step of symmetry restoration recovers also 
the linear properties of the many-body Schrodinger equation. 

The relation between quantum correlations and the two-step method (also called 
the method of hierarchical approximations) is portrayed by the downward vertical arrow 
on the right of figure [l] Indeed, the correlation energy is defined |66] as the difference 
between the restricted Hartree-Fock and exact ground-state energies, i.e., 

-Ecorr = -E'rHF — -^EXD- (1-3) 

As seen from figure [T| starting with the broken-symmetry UHF solution, each further 
approximation captures successively a larger fraction of the correlation energy (1.3); a 



specific example of this process is given in figure [5] below (in section 2.2). 

An alternative approach for studying the emergence of crystalline structures is the 



exact-diagonalizaion method that will be discussed in detail in section 4.1 Like the 



projected wave functions, the EXD many-body wave functions preserve of course all the 
symmetries of the original Hamiltonian. As a result, the intrinsic, or hidden, point-group 
symmetry associated with particle localization and molecule formation is not explicit, 
but it is revealed through inspection of CPDs [one simply uses the exact-diagonalization 
wave function $^™(ri, r2, . . . , r^) in Equation (1.1)] and Equation (1.2), or recognized 
via characteristic trends in the calculated excitation spectra. When feasible, the EXD 
results provide a definitive answer in terms of numerical accuracy, and as such they serve 
as a test to the results obtained through approximation methods (e.g., the above two- 
step method). However, the underlying physics of electron or boson molecule formation 
is less transparent when analyzed with the exact-diagonalization method compared to 
the two-step approach. Indeed, many exact-diagonalization studies of 2D quantum dots 
and trapped bosons in harmonic traps have focused simply on providing high accuracy 
energetics and they omitted calculation of CPDs. However, the importance of using 
CPDs as a tool for probing the many-body wave functions cannot be overstated. For 
example, while exact-diagonalization calculations for bosons in the lowest Landau level 
have been reported rather early [HZl EHl EHl EHl E], the analysis in these studies did 
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not include calculations of the CPDs, and consequently formation of rotating boson 
molecules and particle "crystallization" was not recognized (for further discussion of 
these issues, see Romanovsky et al [601 E3] and Baksmaty et al [72]). 

From the above, it is apparent that both methods, i.e., the two-step method and 
the exact-diagonalization one, complement each other, and it is in this spirit that we 
use them in this report. 

1.6. Experimental signatures of quantum correlations 

Historically, the isolation of a small number [N < 20) of electrons down to a single 
electron was experimentally realized in the so-called "vertical" quantum dots |2]. The 
name vertical QDs derives from the fact that the leads and voltage gates are located in 
a vertical arrangement, on top and below the two-dimensional dot. At zero magnetic 
field, experimental measurements [21 [73] of addition energies, 

A'^En = /iiv+i - fJ^N, (1-4) 

where the chemical potential fiN = En — E^^i, indicated that correlation effects at zero 
and low B are rather weak in such dots, a property that later was attributed to the 
strong screening of the Coulomb interaction in these devices. The measured addition 
energies exhibited maxima at closed electronic shells (A^ = 2, 6, 12, . . .) and at mid-shells 
(A^ = 4, 9, . . .) in agreement with a 2D-harmonic-oscillator central-mean-field model and 
the Hund's rules, and in analogy with the Aufbau principle and the physics of natural 
3D atoms. It was found that the measured ground-state energy spectra for low magnetic 
fields could be understood on the basis of a simple "constant-interaction" model where 
the effect of the two-body Coulomb interaction is reduced phenomenologically to an 
overall classical capacitance, C, characterizing the charging energy Z^e^/(2C) of the 
quantum dot. 

As a result of screening, strong correlation effects and formation of Wigner 
molecules can be expected to occur in vertical dots particularly under the influence 
of high magnetic fields. Evidence about the formation of Wigner molecules in vertical 
quantum dots has been provided recently in Ref. [71], where measured ground-state 
spectra as a function of B for A^ = 3e and N = 4e were reanalyzed with exact- 
diagonalization calculations that included screening. At the time of submission of this 
report, a second ground-state crossing at high B due to strong correlations was also 
demonstrated experimentally in a two-electron vertical quantum dot with an external 
confinement that was smaller than the previously used ones [75] . 

Early theoretical work [20j at zero magnetic field using simply the symmetry broken 
UHF solutions suggested that an unscreened Coulomb repulsion may result in a violation 
of Hund's rules. However, following the two-step method of Refs. [2011211 |22ll23ll2ll25] , 
it has been shown [TB] most recently that the companion step of symmetry restoration 
recovers the Hund's rules in the case of A^ = 4e. 

In addition, the B = results of Ref. [20] suggested that both the maxima of the 
addition energies at closed shells and at mid-shells become gradually weaker (and they 
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eventually disappear) as the strength of the Coulomb interaction (and consequently the 
strength of correlations) increases, leading to formation of "strong" Wigner molecules. 
The qualitative trend of formation of strong Wigner molecules obtained from a relatively 
simple UHF calculation at B = has been confirmed later by more accurate EXD [50|[77] 
and quantum Monte Carlo [TS] calculations, as well as through symmetry restoration 
calculations |23l[Z6], although its experimental demonstration remains still a challenge. 

A more favorable experimental configuration for the development and observation of 
strong interelectron correlations is the so-called "lateral" dot, where the leads and gates 
are located on the sides of the dot and thus screening effects are reduced. Tunability 
of these dots down to a single electron has been achieved only in the last few years 
[79] . Most recently, continually improving experimental techniques have allowed precise 
measurements of excitation spectra of 2e lateral (and anisotropic) quantum dots at zero 
and low magnetic fields [HI HSl EQ]. As discussed in detail in section [sj the behavior of 
these excitation spectra [HI |15] as a function of B provides unambiguous signatures for 
the presence of strong correlations and the formation of Wigner molecules. 

Experimentally observed behavior of two electrons in lateral double QDs [81] 
provides further evidence for strong correlation phenomena. Indeed, instead of 
successively populating delocalized states over both QDs according to a molecular- 
orbital scheme, the two electrons localize on the individual dots according to a Heitler- 
London picture [82]. Theoretically, such strongly correlated phenomena in double 



quantum dots were described in Refs. [201 l22l US]; see section 2.1.4 below. 

Correlations are expected to influence not only the spectral properties of quantum 
dots, but also to effect transport characteristics. Indeed correlation effects may 
underlie the behavior of the transmission amplitudes (magnitude and phase) of an 
electron tunneling through a quantum dot. Such transmission measurements have been 
performed using Aharonov-Bohm interferometry ^3], and an interpretation involving 
strongly correlated states in the form of Wigner molecules has been proposed recently 
[8^ . The quantity that links transport experiments with many-body theory of electrons 
in QDs is the overlap between many-body states with — 1 and A^ electrons, i.e., 
($(A^ — l)|cj|$(A^)), where cj annihilates the jth electron. 

The strength of correlations in quantum dots at zero B can be quantified by the 
Wigner parameter Rw, which is the ratio between the strength of the Coulomb repulsion 
and the one-electron kinetic energy (see section 2.1.2). Naturally, for the case of neutral 
repelling bosons, the corresponding parameter is the ratio between the strength of the 
contact interaction and the one-particle kinetic energy in the harmonic trap, and it 
is denoted as Rs- Larger values of these parameters {Rw or Rg) result in stronger 
correlation effects. 

Progress in the ability to experimentally control the above parameters has been 
particularly impressive in the case of ultracold trapped bosons. Indeed, realizations of 
continuous tunability of Rs over two orders of magnitude (from 1 to 5 [85] and from 
5 to 200 |Sn]) has been most recently reported in quasi-linear harmonic traps. Such 
high values of Rs allowed experimental realization of novel strongly correlated states 
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drastically different from a Bose-Einstein condensate. This range of high values of Rs is 
known as the Tonks- Girardeau regime and the corresponding states are one-dimensional 
analogues of molecular structures made out of localized bosons. In two dimenional traps, 
it has been predicted that such large values of Rs lead to the emergence of crystalhne 
phases [201 E3] • 

The high experimental control of optical lattices has also been exploited for the 
creation [87J of novel phases of ultracold bosons analogous to Mott insulators; such 
phases are related to the formation of electron puddles discussed in section 2.1.4 and to 
the fragmentation of Bose-Einstein condensates [88] . 



1.7. Plan of the report 

The plan of the report can be visualized through the table of contents. Special attention 
has been given to the Introduction, which offers a general presentation of the subject 
of symmetry breaking and quantum correlations in confined geometries - including a 
discussion of the differences with the case of extended systems, a historical background 
from other fields, and a diagrammatic synopsis of the two-step method of symmetry 
breaking/symmetry restoration. 

The theoretical framework and other technical methodological background are 
presented in Section 2 (symmetry breaking/symmetry restoration in quantum dots). 
Section 3 (symmetry breaking/symmetry restoration for trapped ultracold bosons), and 
Section 4 (exact-diagonalization approaches). Section 4 includes also a commentary on 
quantum Monte Carlo methods. 

For the case of semiconductor quantum dots, the main results and description of 
the strongly correlated regime are presented in Sections 5, 6, and 7, with Section 5 
focusing on the case of two electrons and its historical significance. Section 8 is devoted 
to a description of the strongly-correlated regime of trapped repelling bosons. 

Finally, a summary is given in Section 9, and the Appendix offers an outline of the 
Darwin-Fock single-particle spectra for a two-dimensional isotropic oscillator under a 
perpendicular magnetic field or under rotation. 

We note that the sections on trapped bosons (Section 3 and Section 8) can be read 
independently from the sections on quantum dots. 



2. Symmetry breaking and subsequent symmetry restoration for electrons 
in confined geometries: Theoretical framework 

The many-body Hamiltonian describing N electrons confined in a two-dimensional QD 
and interacting via a Coulomb repulsion is written as 

N N N 2 

^ = E^(0 + EE— • (2-1) 

1=1 1=1 j>i ^'"iJ 



In Equation (2.1), k is the dielectric constant of the semiconducting material and 
fij = ki — rj|. The single-particle Hamiltonian in a perpendicular external magnetic 
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field B is given by 

where the external confinement is denoted by V{x,y), the vector potential A is given 
in the symmetric gauge by 

A(r) = is X r = ^(-Sy, Bx, 0), (2.3) 



and the last term in (2.2) is the Zeeman interaction with g* being the effective Lande 
factor, fiB the Bohr magneton, s the spin of an individual electron and m* is the 
effective electron mass. The external potential confinement V{x, y) can assume various 
parametrizations in order to model a single circular or elliptic quantum dot, or a 
quantum dot molecule. Of course, in the case of an elliptic QD, one has 

V{x,y) = ^m*{ujy + ujy), (2.4) 

which reduces to the circular QD potential when = ujy = ujq. The appropriate 
parametrization of V{x,y) in the case of a double QD is more complicated. In our 
work, we use a parametrization based on a 2D two-center oscillator with a smooth 
necking. This latter parametrization is described in detail in Refs. [23t H6] . where 
readers are directed for further details. In contrast with other parametrizations based 
on two displaced inverted Gaussians |89], the advantage of the two-center oscillator is 
that the height of the interdot barrier, the distance between the dots, the ellipticity of 
each dot, and the gate potentials of the two dots (i.e., the relative potential wells in the 
neighboring dots) can be varied independently of each other. 

A prefactor multiplying the Coulomb term in Equation (2.1) (being either an 
overall constant 7 as in section bA_ below, or having an appropriate position-dependent 
functional form |321I1S]) is used to account for the reduction of the Coulomb interaction 
due to the finite thickness of the electron layer and to additional screening (beyond that 
produced by the dielectric constant of the material) arising from the formation of image 
charges in the gate electrodes |90] . 



2.1. Mean- field description and unrestricted Hartree-Fock 

Vast literature is available concerning mean-field studies of electrons in quantum 
dots. Such publications are divided mainly into applications of density functional 
theory P EH El El El ESI ES] and the use of Hartree-Fock methods [201 EZl 
ESI EH ES Ea [rnni [MI [M [M HDI]. The latter include treatments according to 
the restricted Hartree-Fock [97], unrestricted Hartree-Fock with spin, but not space, 
symmetry breaking [HSl ESI HOOj - unrestricted Hartree-Fock with spin and/or space 
symmetry breaking |20l |25l [271 El llOlj I102[ 1103] , and the so-called Brueckner Hartree- 
Fock [TMfTUS] . 

From the several Hartree-Fock variants mentioned above, only the UHF with 
consideration of both spin and space symmetry unrestrictions has been able to describe 
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formation of Wigner molecules, and in the following we will exclusively use this 
unrestricted version of Hartree-Fock theory. The inadequacy of the density-functional 



theory in describing Wigner molecules will be discussed in section 2.3 



2.1.1. The self- consistent Pople-Nesbet equations. The unrestricted Hartree-Fock 
equations used by us are an adaptation of the Pople-Nesbet |106j equations described 
in detail in Ch 3.8 of Ref. [19]. For completeness, we present here a brief description of 
these equations, along with pertinent details of their computational implementation by 
us to the 2D case of semiconductor QDs. 

We start by requesting that the unrestricted Hartree-Fock many-body wave function 
for electrons is represented by a single Slater determinant 

^uhf(xi, . . . , xtv) = -^det[xi(xi), X2(x2), • • • , Xn{^n], (2.5) 

where ^'^'^ a set of spin orbitals, with the index x denoting both the space 

and spin coordinates. Furthermore, we take Xi(x) = ipi{r)a for a spin-up electron 
and = foi' ^ spin-down electron. As a result, the UHF determinants 

in this report are eigenstates of the projection of the total spin with eigenvalue 
5*2 = (A^" — N^)/2, where A^"(^) denotes the number of spin up (down) electrons. 
However, these Slater determinants are not eigenstates of the square of the total spin, 
S^, except in the fully spin polarized case. 

According to the variational principle, the best spin orbitals must minimize the total 
energy (\I/uhf|'^|^uhf)- By varying the spin orbitals [Xj(x)] under the constraint that 
they remain orthonormal, one can derive the UHF Pople-Nesbet equations described 
below. 

A key point is that electrons with a (up) spin will be described by one set of spatial 
orbitals {ipflJ = 1,2, . . . , K}, while electrons with (3 (down) spin are described by a 
different set of spatial orbitals = 1,2, ... , K}; of course in the restricted Hartree- 

Fock = = ^j. Next, one introduces a set of basis functions {ffilfJ^ = 1, 2, ... , K} 
(constructed to be orthonormal in our 2D case), and expands the UHF orbitals as 

C = EC^>., ^ = l,2,...,K, (2.6) 

^f = EC^5^^, ^ = l,2,...,K. (2.7) 
11=1 

The UHF equations are a system of two coupled matrix eigenvalue problems 
resolved according to up and down spins, 

F"^C" = C"E" (2.8) 

p/3a(^/3 ^ c^-E^, (2.9) 

where F"'^^^"-' are the Fock-operator matrices and C"'-^-' are the vectors formed with the 



coefficients in the expansions (2.6) and (2.7). The matrices E"*^^) are diagonal, and as 
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a result equations (2.8) and (2.9) are canonical (standard). Notice that noncanonical 
forms of HF equations are also possible (see Cli 3.2.2 of Ref. |19j|). Since the self- 
consistent iterative solution of the HF equations can be computationally implemented 
only in their canonical form, canonical orbitals and solutions will always be implied, 
unless otherwise noted explicitly. We note that the coupling between the two UHF 



equations (2.8) and (2.9) is given explicitly in the expressions for the elements of the 



Fock matrices below [(2.12) and (2.13)]. 

Introducing the density matrices P"(^) for a{P) electrons. 



Nf3 



where A^" + A^'' = A^, the elements of the Fock-operator matrices are given by 
= i/,. + EE^r.[(/i^kA) - (/io-IAz.)] +EE^il(/^^I^A) 



(2.10) 



(2.11) 



(2.12) 



F^^ = H,, + J2J2P^^[i^^a\uX) - (/xa|Az/)] + E E ^^.(/^^kA), (2.13) 



fii/ 



where H^^, are the elements of the single electron Hamiltonian (with an external 
magnetic field B and an appropriate potential confinement), and the Coulomb repulsion 
is expressed via the two-electron integrals 

dridr2</?*(ri)(/?;(r2) 



(^(t|z/A) 



e 



(2.14) 



|ri - r2| 

with K being the dielectric constant of the semiconductor material. Of course, the Greek 
indices fi, z/. A, and a run from 1 to K. 



The system of the two coupled UHF matrix equations (2.8) and (2.9) is solved 



self consistently through iteration cycles. For obtaining the numerical solutions, we have 
used a set of K basis states ipi's that are chosen to be the product wave functions formed 
from the eigenstates of one-center (single QD) and/or two-center [221 SS] (double QD) 
one-dimensional oscillators along the x and y axes. Note that for a circular QD a value 
K = 78 corresponds to all the states of the associated 2D harmonic oscillator up to and 
including the 12th major shell. 

The UHF equations preserve at each iteration step the symmetries of the many-body 
Hamiltonian, if these symmetries happen to be present in the input (initial) electron 
density of the iteration (see section 5.5 of Ref. [IS]). The input densities into the 
iteration cycle are controlled by the values of the P"^ and Pf^ matrix elements. Two 
cases arise in practice: (i) Symmetry adapted RHF solutions are extracted from (2.8) 



and (2.9) by using as input P^^ = P{;=0 for the case of closed shells (with or without 
an infinitesimally small B value). For open shells, one needs to use an infinitesimally 
small value of B. With these choices, the output of the first iteration (for either closed 
or open shells) is the single-particle spectrum and corresponding electron densities at 
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B = associated with the Hamiltonian in ( 2.2[ ) (the small value of B mentioned above 



guarantees that the single-particle total and orbital densities are circular), (ii) For 
obtaining broken-symmetry UHF solutions, the input densities must be different in an 
essential way from the ones mentioned above. We have found that the choice P"^ = 1 and 
Pf^ = usually produces broken-symmetry solutions (in the regime where symmetry 
breaking occurs). 

Having obtained the selfconsistent solution, the total UHF energy is calculated as 

^UHF = I E E[(^.; + K)H,u + P:,F^J + P'.X^]. (2.15) 

We note that the Pople-Nesbet UHF equations are primarily employed in Quantum 
Chemistry for studying the ground states of open-shell molecules and atoms. Unlike our 
studies of QDs, however, such chemical UHF studies consider mainly the breaking of 
the total spin symmetry, and not that of the space symmetries. As a result, for purposes 
of emphasis and clarity, we have often used (see, e.g., our previous papers) prefixes to 
indicate the specific unrestrictions (that is removal of symmetry restrictions) involved 
in our UHF solutions, i.e., the prefix s- for the total-spin and the prefix S- for the space 
unrestriction. 

The emergence of broken-symmetry solutions is associated with instabilities of the 
restricted HF solutions, i.e., the restricted HF energy is an extremum whose nature as a 
minimum or maximum depends on the positive or negative value of the second derivative 
of the HF energy. The importance of this instability problem was first highlighted in a 
paper by Overhauser [107] . Soon afterwards, the general conditions for the appearance of 
such instabilities (analyzed within linear response and the random-phase approximation) 
were discussed by Thouless in the context of nuclear physics |108] . Subsequently, 
the Hartree-Fock stability/instability conditions were re-examined |109l 1110] . using a 
language from (and applications to) the field of quantum chemistry. For comprehensive 
reviews of mean-field symmetry breaking and the Hartree-Fock methods and instabilities 
in the context of quantum chemistry, see the collection of papers in Ref . jlll] . 



2.1.2. The Wigner parameter and classes of spontaneous symmetry breaking solutions. 
Using the self-consistent (spin-and-space) unrestricted Hartree-Fock equations presented 
in the previous section, we found [20], for zero and low magnetic fields, three classes 
of spontaneous symmetry breakings in circular single QDs and in lateral quantum dot 
molecules (i.e., formation of ground states of lower symmetry than that of the confining 
potentials). These include the following: 

(I) Wigner molecules in both QDs and quantum dot molecules, i.e., (spatial) 
localization of individual electrons within a single QD or within each QD comprising 
the quantum dot molecule. 

(II) formation of electron puddles in quantum dot molecules, that is, localization 
of the electrons on each of the individual dots comprising the quantum dot molecule, 
but without localization within each dot, and 
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Figure 2. (Color online) UHF electron density in a parabolic QD for TV = 19 and 
Sz = 19/2, exhibiting breaking of the circular symmetry at Rw = 5 and B = 0. The 
choice of the remaining parameters is: hujo = 5 meV and m* = 0.067me. Distances 
are in nanometers and the electron density in 10^* nm^^. 

(Ill) pure spin-density waves (SDWs) which are not accompanied by spatial 
locahzation of the electrons j91j . 

It can be shown that a central-mean-field description (associated with the RHF) at 
zero and low magnetic fields may apply in the case of a circular QD only for low values 
of the Wigner parameter 

/?w = Q/huo, (2.16) 

where Q is the Coulomb interaction strength and huJo is the energy quantum of the 
harmonic potential confinement (being proportional to the one-particle kinetic energy); 
Q = e^/(K/o), with k being the dielectric constant, /q = {h/{m*uJo)Y^'^ ^^e spatial 
extension of the lowest state's wave function in the harmonic (parabolic) confinement, 
and m* the effective electron mass. 

Furthermore, we find that Wigner molecules (SSB class I) occur in both QDs and 
quantum dot molecules for i?w > 1- Depending on the value of Rw, one may distinguish 
between "weak" (for smaller Rw values) and "strong" (for larger Rw values) Wigner 
molecules, with the latter termed sometimes as "Wigner crystallites" or "electron 
crystallites." The appearance of such crystalline structures may be regarded as a 
quantum phase transition of the electron liquid upon increase of the parameter Rw- 
Of course, due to the finite size of QDs, this phase transition is not abrupt, but it 
develops gradually as the parameter Rw varies. 

For quantum dot molecules with i?w < 1, Wigner molecules do not develop and 
instead electron puddles may form (SSB class II). For single QDs with i?w < 1, we find 
in the majority of cases that the ground-states exhibit a central-mean-field behavior 
without symmetry breaking; however, at several instances (see an example below), a 
pure SDW (SSB class III) may develop. 
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Figure 3. (Color online) UHF ground-state of a 6e quantum dot molecule (double 
dot), with parameters resulting in formation of two non-crystallized electron puddles 
(akin to dissociation of the quantum dot molecule in two QDs with 3 electrons each). 
Left: total electronic density. Right: contour plots of the densities (orbital squares) 
of the three individual orbitals localized on the left dot, with spin polarization of 
the orbitals as indicated. The choice of parameters is: Tiloq = 5 meV (harmonic 
confinement of each dot), d = 70 nm (distance berween dots), Vi, = 10 meV (interdot 
barrier), m* = 0.067me (electron effective mass), and k = 20 (dielectric constant). 
Lengths {x and y axes) in nm, density distribution (vertical axis) in 10^^ nm~^. 



2.1.3. Unrestricted Hartree-Fock solutions representing Wigner molecules. As a typical 
example of a Wigner-molecule solution that can be extracted from the UHF equations, 
we mention the case of = 19 electrons for hujQ = 5 meV, Rw = 5 (k = 3.8191), and 
B = 0. Figure [2] displays the total electron density of the broken-symmetry UHF 
solution for these parameters, which exhibits breaking of the rotational symmetry. 
In accordance with electron densities for smaller dot sizes published by us earlier 
[20l |2T] the electron density in figure |2] is highly suggestive of the formation of a 
Wigner molecule, with a (1,6,12) ring structure in the present case; the notation 
{ni,n2, . . . ,nr) signifies the number of electrons in each ring: ni in the first, n2 in 
the second, and so on. This polygonal ring structure agrees with the classical one 
(that is the most stable arrangement of 19 point charges in a 2D circular harmonic 
confinement \112\ \113\ I114| j|]). and it is sufficiently complex to instill confidence that 
the Wigner-molecule interpretation is valid. The following question, however, arises 
naturally at this point: is such molecular interpretation limited to the intuition provided 
by the landscapes of the total electron densities, or are there deeper analogies with 
the electronic structure of natural 3D molecules? The answer to the second part 
of this question is in the affirmative. Indeed, it was found [25] that SSB results in 
the replacement of a higher symmetry by a lower one. As a result, the molecular 
UHF solutions exhibit point-group spatial symmetries that are amenable to a group- 
theoretical analysis in analogy with the case of 3D natural molecules. 

II These references presented extensive studies pertaining to the geometrical arrangements of classical 
point charges in a two-dimensional harmonic confinement. 
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2.1.4- Unrestricted Hartree-Fock solutions representing electron puddles. An example 
of formation of electron puddles in quantum dot molecules, that is, localization of the 
electrons on each of the individual dots comprising the quantum dot molecule, but 
without localization within each dot, is presented in figure |3j We consider the case of 
= 6 electrons in a double dot under field-free conditions (5 = 0); with parameters 
hujQ = 5 meV (harmonic confinement of each dot), d = 70 nm (distance berween dots), 
Vb = 10 meV (interdot barrier) and m* = 0.067me (electron effective mass). Reducing 
the Rw value (with reference to each constituent QD) to 0.95 (i.e., for a dielectric 
constant n = 20) guarantees that the ground-state of the 6e quantum dot molecule 
consists of electron puddles [SSB of type II, figure |3]. In this case, each of the electron 
puddles (on the left and right dots) is spin-polarized with total spin projection Sz = 1/2 
on the left QD and Sz = —1/2 on the right QD. As a result, the singlet and triplet states 
of the whole quantum dot molecule are essentially degenerate. Note that the orbitals 
on the left and right dots [see, e.g., those on the left dot in figure |3] (right)] are those 
expected from a central-mean-field treatment of each individual QD, but with slight 
(elliptical) distortions due to the interdot interaction and the Jahn- Teller distortion 
associated with an open shell of three electrons (in a circular harmonic confinement). 
Note the sharp contrast between these central-mean-field orbitals and corresponding 
electron density (figure Isl) with the electron density and the three orbitals associated 



with formation of a Wigner molecule inside a single QD [see, e.g., figure[6jin section 2.2.2 
below] . 

The formation of electron puddles described above can be also seen as a form of 
dissociation of the quantum dot molecule. We found that only for much lower values 
of Rw (< 0.20, i.e., k > 90.0) the electron orbitals do extend over both the left and 
right QDs, as is usually the case with 3D natural molecules (molecular-orbital theory). 
Further examples and details of these two regimes (dissociation versus molecular-orbital 
description) can be found in Refs. [221 SH] • 



2.1.5. Unrestricted Hartree-Fock solutions representing pure spin density waves within 
a single quantum dot. Another class of broken-symmetry solutions that can appear in 
single QDs are the spin density waves. The SDWs are unrelated to electron localization 
and thus are quite distinct from the Wigner molecules [20]; in single QDs, they 
were obtained [H] earlier within the framework of spin density functional theory. To 
emphasize the different nature of spin density waves and Wigner molecules, we present 
in figure |4] an example of a SDW obtained with the UHF approach [the corresponding 
parameters are: A^ = 14, Sz = 0, i?w = 0.8 (k = 23.8693), and B = 0]. Unlike the case 
of Wigner molecules, the SDW exhibits a circular electron density [see figure |4]^ a)], and 
thus it does not break the rotational symmetry. Naturally, in keeping with its name, 
the SDW breaks the total spin symmetry and exhibits azimuthal modulations in the 
spin density [see figure |4](b); however, the number of humps is smaller than the number 
of electrons] . 

We mention here that the possibility of ground-state configurations with uniform 
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Figure 4. (Color online) UHF solution in a parabolic QD exhibiting a pure spin 
density wave for N = 14, Sz = 0, i?w = 0.8, and B — 0. (a) The total electron density 
exhibiting circular symmetry; (b) The spin density exhibiting azimuthal modulation 
(note the 12 humps whose number is smaller than the number of electrons; on the 
contrary in the case of a Wigner molecule, the number of humps in the electron density 
is always equal to N). The choice of the remaining parameters is: Kluq = 5 meV and 
m* = 0.067me. Distances are in nanometers and the electron (ED) and spin (SD) 
densities in 10^* nm~^. 



electron density, but nonuniform spin density, was first discussed for 3D bulk metals 
using the HF method in Ref. |115] . 

The SDWs in single QDs appear for i?w < 1 and are of lesser importance; thus in the 
following we will exclusively study the case of Wigner molecules. However, for ^ 1; 



formation of a special class of SDWs (often called electron puddles, see section 2.1.4) 
plays an important role in the coupling and dissociation of quantum dot molecules (see 
Ref. and Ref. 

2.2. Projection techniques and post-Hartree-Fock restoration of broken symmetries 



As discussed in section 1.5, for finite systems the symmetry broken UHF solutions are 
only an intermediate approximation. A subsequent step of post-Hartree-Fock symmetry 
restoration is needed. Here we present the essentials of symmetry restoration while 
considering for simplicity the case of two electrons in a circular parabolic QD. 

Results obtained for various approximation levels for a two-electron QD with B = 
and Rw = 2.40 (that is, in the Wigner-molecule regime) are displayed in figure |5} In 
these calculations [23], the spin projection was performed following reference [31], i.e., 
one constructs the wave function 

^Spin-p(s) = Pspin(s)^UHF, (2.17) 

where ^E'uhf is the original symmetry-broken UHF determinant (which is already by 



construction an eigenstate of the projection Sz of the total spin). In ( |2.17 ), the spin 



projection operator (projecting into a state which is an eigenstate of the square of the 
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Figure 5. (Color online) Various approximation levels for the lowest singlet state 
of a field- free two-electron QD with R\y = 2 AO. The corresponding energies (in 
meV) are shown at the bottom of the figure, (a): Electron density of the RHF 
solution, exhibiting circular symmetry (due to the imposed symmetry restriction). The 
correlation energy Scorr = 2.94 meV, is defined as the difference between the energy of 
this state and the exact solution [shown in frame (e)]. (bl) and (b2): The two occupied 
orbitals (modulus square) of the symmetry-broken "singlet" UHF solution (bl), with 
the corresponding total electron density exhibiting non-circular shape (b2) . The energy 
of the UHF solution shows a gain of 44.3% of the correlation energy, (c): Electron 
density of the spin-projected singlet (Spin-P), showing broken spatial symmetry, but 
with an additional gain of correlation energy, (d): the spin-and-angular-momentum 
projected state (S&AMP) exhibiting restored circular symmetry with a 73.1% gain 
of the correlation energy. The choice of parameters is: dielectric constant k = 8, 
parabolic confinement hujg = 5 meV, and effective mass m* = 0.067me. Distances are 
in nanometers and the densities in 10^^ nm^^. 



total spin) is given by 



n 



S2 - S'{S' + 1) 



(2.18) 



where the index s' runs over the quantum numbers associated with the eigenvalues 
s'{s' + 1) of (in units of h^), with S being the total spin operator. For two electrons, 



the projection operator reduces to V, 



s,t 
spin 



1 T^i2, where the operator interchanges 



the spins of the two electrons; the upper (minus) sign corresponds to the singlet (s 
supersript), and the lower (plus) sign corresponds to the triplet (t superscript) state. 

The angular momentum projector (projecting into a state with total angular 
momentum L) is given by 

27rVL = /'"d7exp[-i7(L - L)], (2.19) 
Jo 

where L = li + 12 is the total angular momentum operator. As seen from 
(2.19), apphcation of the projection operator Vl to the spin-restored state ^'spin-p('S) 



corresponds to a continuous configuration interaction expansion of the wave function 



that uses, however, non- orthogonal orbitals (compare section 4.1) 



The application of the projection operator Vl to the state \E'spin-p(s) generates a 
whole rotational band of states with good angular momenta (yrast band). The energy 
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of the projected state with total angular momentum L is given by 

EpRj(L)=/ /i(7)e'^^d7/ / n(7)ei^M7, (2.20) 



/o / 

with /i(7) = (^'spin-p(s;0)|7^|^'spin-p(s;7)) and n(7) = (^spm-p(s; 0)|^spin-p(s; 7)), 
where \l/spin-p('S; 7) is the spin-restored (i.e., spin-projected) wave function rotated by 
an azimuthal angle 7 and H is the many-body Hamiltonian. We note that the UHF 
energies are simply given by -Euhf = /i(0)/'n.(0). 

In the following we focus on the ground state of the two-electron system, i.e., 
L = 0. The electron densities corresponding to the initial RHF approximation [shown 
in figure [S^a)] and the final spin-and-angular-momentum projection (S&AMP) [shown in 
figure [sj^d)], are circularly symmetric, while those corresponding to the two intermediate 
approximations, i.e., the UHF and spin-projected solutions [figure |5](b 2) and figure [5|c), 
respectively] break the circular symmetry. This behavior illustrates graphically the 
meaning of the term "restoration of symmetry," and the interpretation that the UHF 
broken-symmetry solution refers to the intrinsic (rotating) frame of reference of the 
electron molecule. In light of this discussion the final projected state is called a rotating 
electron or (Wigner) molecule. 



Expressions (2.19) and (2.20) apply directly to REMs having a single polygonal ring 
of localized electrons, with L = J2^=i h- For a generalization to electron molecules 



with multiple concentric polygonal rings, see section |2.2.1| below. 

For restoring the total spin, an alternative method to the projection formula (2.18) 
can be found in the literature [33]. We do not make use of this alternative formulation 
in this report, but we briefly describe it here for the sake of completeness. Based on the 
formal similarity between the 3D angular momentum and the total spin, one can apply 
the formula by Peierls and Yoccoz [30] and obtain the projection operator 

2s + 1 



^1.. = J dTDg^{T)n{T), (2.21) 

where Dg'^glT) are the 3D Wigner D functions |116j . F is a shorthand notation for the 
set of the three Euler angles {(f),9,ilj), and 

7^(^) = e-i^s^e-iosy^-i^s. (2.22) 



is the rotation operator in spin space. In (2.21), the indices of the Wigner D functions 
are s, Sz, and q. 

The operator Vg^^ extracts from the symmetry broken wave function a state with 
a total spin S and projection Sz along the laboratory z axis. However, q is not a 
good quantum number of the many-body Hamiltonian, and the most general symmetry 
restored state is written as a superposition over the components of q, i.e., 

^Spin-P(s, Sz-, l) = T.aiVl^'^VKF, (2.23) 

g 

where the coefficients are determined through a diagonalization of the many-body 
Hamiltonian in the space spanned by the nonorthogonal PI^q^I/uhf (see also Refs. 
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[1171 IllSj ). In (2.23), the index i reflects the possible degeneracies of spin functions 
with a given good total-spin quantum number s |119j . which is not captured by (2.18). 

The Peierls-Yoccoz formulation for recovering spin-corrected wave functions applies 
also in the case when the UHF determinants violate in addition the conservation of spin 
projection |33], unlike the projector Pspml-s) [see (2.18)] which acts on UHF determinants 
having a good 5*2 = (A^" — N^)/2 according to the Pople-Nesbet theory presented in 
section I2.1.1[ 

In the literature [18], there are two distinguishable implementations of symmetry 
restoration: variation before projection (VBP) and variation after projection (YAP). In 
the former, which is the one that we mostly use this report, mean-field solutions with 
broken symmetry are first constructed and then the symmetry is restored via projection 
techniques as described above. In the latter, the projected wave function is used as the 
trial wave function directly in the variational principle (in other words the trial function 
is assured to have the proper symmetry). 

The VAP is in general more accurate, but more difficult to implement numerically, 
and it has been used less often in the nuclear-physics literature. In quantum chemistry, 
the generalized valence bond method [12U[ . or the spin-coupled valence bond method 
[121] . describing covalent bonding between pairs of electrons, employ a variation after 
projection. 

For quantum dots, the variation after projection looks promising for reducing the 
error of the VBP techniques in the transition region from mean-field to Wigner-molecule 
behavior, where this error is the largest. In fact, it has been found that the discrepancy 
between variation-before-projection techniques and exact solutions is systematically 
reduced [231 ESI [122] for stronger symmetry breaking (increasing Ry/ and/or increasing 
magnetic field). 

Moreover, in the case of an applied magnetic field (quantum dots) or a rotating trap 
(Rose gases), our VBP implementation corresponds to projecting cranked symmetry- 
unrestricted Slater determinants [123] . This is because of the "cranking" terms 
—'hujcL/2 or —hQL that contribute to the many-body Hamiltonian Ti, respectively, 
with uJc = eB/{m*c) being the cyclotron frequency and Vt the rotational frequency of 
the trap; these terms arise in the single-particle component of Ti [see Equation (2.2) 
in section [2] and Equation (8.3) in section [s]. The cranking form of the many-body 
Hamiltonian is particularly advantageous to the variation before projection, since the 
cranking method provides a first-order approximation to the variation- a/ier-projection 
restoration of the total angular- momentum L [124] (see also Ch 11.4.4 in Ref. [T8]). 



2.2.1. The REM microscopic method in medium and high magnetic field. In our 
method of hierarchical approximations, we begin with a static electron molecule, 
described by an unrestricted Hartree-Fock determinant that violates the circular 
symmetry [201 1231 12S] • Subsequently, the rotation of the electron molecule is described by 
a post-Hartree-Fock step of restoration of the broken circular symmetry via projection 
techniques [221 1231 1211 1251 [261 [511 [53]. Since we focus here on the case of strong i?. 
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we can approximate the UHF orbitals (first step of our procedure) by (parameter free) 
displaced Gaussian functions; that is, for an electron localized at Hj (Zj), we use the 
orbital l53] 



1 / \z-Z, 



|2 



uiz,Z,) = — exp -^—>^-i^{z,ZfB) , (2.24) 



with A = / = ^h/m*uj; uj = ^uJq + u'^/A, where oOc = eB/{m*c) is the cyclotron 
frequency and uo specifies the external parabolic confinement. We have used complex 
numbers to represent the position variables, so that z = x + iy, Zj = Xj + iYj. The 
phase guarantees gauge invariance in the presence of a perpendicular magnetic field and 
is given in the symmetric gauge by ip{z, Zj; B) = {xYj — yXj)/2l'^, with Ib = ^Jhc/eB. 

For an extended 2D system, the Zj^s form a triangular lattice [SHI I125j . For finite 
N, however, the Z/s coincide [211 |26l EU |52l |53] with the equilibrium positions [forming 
r concentric regular polygons denoted as (ni, n2, . . . , n^)] of = Ylg=i classical point 
charges inside an external parabolic confinement |114] . In this notation, ni corresponds 
to the innermost ring with ni > 0. For the case of a single polygonal ring, the notation 
(0, N) is often used; then it is to be understood that ni = N. 

The wave function of the static electron molecule is a single Slater determinant 
|v|>UHF|-^j^ made out of the single-electron wave functions u{zi,Zi), i = 1,...,N. 
Correlated many-body states with good total angular momenta L can be extracted 
[211 I2ni EH E3] (second step) from the UHF determinant using projection operators. 
The projected rotating electron molecule state is given by 

1$!™) = r... r • • • d7.|^"""(7i, • • • , 7.)) exp (i ± 7,L, 

Jo Jo y 

(2.25) 

Here L = J2g^i Lq and |\E''^^^[7]) is the original Slater determinant with all the single- 
electron wave functions of the qth ring rotated (collectively, i.e., coherently) by the same 



azimuthal angle 7g. Note that (2.25) can be written as a product of projection operators 
acting on the original Slater determinant [i.e., on |\E'^^^(7i = 0, . . . ,7,. = 0))]. Setting 
A = IbV^ restricts the single-electron wave function in (2.24) to be entirely in the lowest 
Landau level (see Appendix in Ref. [53]). The continuous-configuration- interaction 
form of the projected wave functions [i.e., the linear superposition of determimants in 



(2.25)] implies a highly entangled state. We require here that B is sufficiently strong so 
that all the electrons are spin-polarized and that the ground-state angular momentum 
L > Lq = J2i^^ i = N(N — l)/2 (or equivalently that the fractional filling factor 
u = Lq/L < 1). The state corresponding to Lq is a single Slater determinant in the 
lowest Landau level and is called the "maximum density droplet" |126j . For high B, the 
calculations in this paper do not include the Zeeman contribution, which, however, can 
easily be added (for a fully polarized dot, the Zeeman contribution to the total energy 
is Ng* fiBB /2, with g* being the effective Lande factor and ^b the Bohr magneton). 

Due to the point-group symmetries of each polygonal ring of electrons in the UHF 
wave function, the total angular momenta L of the rotating crystalline electron molecule 
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are restricted to the so-called magic angular momenta, i.e. 

r 

Lm = Lq + kqUq, 



(2.26) 



q=l 



where the /cg's are non-negative integers (when ni = 1, ki = 0). 

Magic angular momenta associated with multiple rings have been discussed in Refs. 
[2T, '261 EH Ea [53]. For the simpler cases of (0, A^) or (1, - 1) rings, see, e.g., Ref. 
[127j and Ref. [HH]. 

The partial angular momenta associated with the gth ring, Lq [see (2.25)], are given 

by 



^0,q + kqTlq, 



(2.27) 



where Lo,q = EiLT,+i(^ - 1) with ig = ELi tt-s {h = 0), and Lq = J2l=i ^o,r 
The energy of the REM state ( |2^ is given [SI EH |52l [53] by 



^REM 



2n 



/i([7])e'W-[i]d[^] 



2tt 



n([7])e^W-Wd[^] 



(2.28) 



with the Hamiltonian and overlap matrix elements ^^([7]) = ([0])|7i|\E' ([7])) 

and n([7]) = (*uhfqq])|^uhfq^])^^ respectively, and [7] • [L] = TJq=ilqLq. The UHF 
energies are simply given by -Euhf = ^([0])/'^([0]). 

The crystalline polygonal-ring arrangement (ni, n2, . . . , n^) of classical point charges 
is portrayed directly in the electron density of the broken-symmetry UHF, since the 
latter consists of humps centered at the localization sites Zj's (one hump for each 
electron). In contrast, the REM has good angular momentum and thus its electron 
density is circularly uniform. To probe the crystalline character of the REM, we use 
the conditional probability distribution (CPD) defined in (1.1). P(r,ro) is proportional 
to the conditional probability of finding an electron at r, given that another electron is 
assumed at Fq. This procedure subtracts the collective rotation of the electron molecule 
in the laboratory frame of referenece, and, as a result, the CPDs reveal the structure of 
the many body state in the intrinsic (rotating) reference frame. 



2.2.2. Group structure and sequences of magic angular momenta. It has been 
demonstrated that the broken-symmetry UHF determinants and orbitals describe 
2D electronic molecular stuctures (Wigner molecules) in close analogy with the case 
of natural 3D molecules. However, the study of Wigner molecules at the UHF level 
restricts their description to the intrinsic (nonrotating) frame of reference. Motivated 
by the case of natural atoms, one can take a subsequent step and address the properties 
of collectively rotating Wigner molecules in the laboratory frame of reference. As is well 
known, for natural atoms, this step is achieved by writing the total wave function of 
the molecule as the product of the electronic and ionic partial wave functions. In the 
case of the purely electronic Wigner molecules, however, such a product wave function 
requires the assumption of complete decoupling between intrinsic and collective degrees 
of freedom, an assumption that might be justifiable in limiting cases only. The simple 



Symmetry breaking and quantum correlations 



28 




Figure 6. (Color online) The UHF solution exhibiting breaking of the circular 
symmetry for iV 3 and Sz — 1/2 at Rw = 10 and B = 0. (a-b): real orbitals 
for the two spin-up electrons, (c): real orbital for the single spin-down electron, (d): 
total electron density (ED), (e): spin density (SD, difference of the spin-up minus 
the spin-down partial electron densities). The choice of the remaining parameters is: 
hajQ = 5 meV and m* — O.OGTmg. Distances are in nanometers. The real orbitals are 
in 10""^ nm"""^ and the densities (electron density and spin density) in 10^* nm"^. The 
arrows indicate the spin direction. 



product wave function was used in earlier treatments of Wigner molecules; see, e.g., Ref. 
|128] . The projected wave functions employed here are integrals over such product wave 
functions, and thus they account for quantal fluctuations in the rotational degrees of 
freedom. The reduction of the projected wave functions to the limiting case of a single 
product wave function is discussed in Ch 11.4.6.1 of Ref. [18] . 

As was discussed earlier, in the framework of the broken-symmetry UHF solutions, 
a further step is needed - and this companion step can be performed by using the 
post-Hartree-Fock method of restoration of broken symmetries via projection techniques 
(see section 2.2). In this section, we use this approach to illustrate through a couple 
of concrete examples how certain universal properties of the exact solutions, i.e., the 
appearance of magic angular momenta |127l I128[ 11291 11301 I13H \132\ I133j in the 
exact rotational spectra, relate to the symmetry broken UHF solutions. Indeed, we 
demonstrate that the magic angular momenta are a direct consequence of the symmetry 
breaking at the UHF level and that they are determined fully by the molecular symmetries 
of the UHF determinant. 

As an illustrative example, we have chosen the relatively simple, but non trivial 
case, of iV = 3 electrons. For B = 0, both the = 1/2 and Sz = 3/2 polarizations 
can be considered. We start with the Sz = 1/2 polarization, whose broken-symmetry 
UHF solution |25] is portayed in figure [6] and which exhibits a breaking of the total 
spin symmetry in addition to the rotational symmetry. Let us denote the corresponding 
UHF determinant [made out of the three spin orbitals in figure [6ta), figure |6](b), and 
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figure [6]^c)] as | ITT)- We first proceed with the restoration of the total spin by noticing 
that I ITT) has a lower point-group symmetry (see Ref. |25]) than the C^v symmetry of 
an equilateral triangle. The C^^ symmetry, however, can be readily restored by applying 



the projection operator (2.19) to | |t|) and by using the character table of the cyclic 
C3 group (see Table I in Ref. j25j). Then for the intrinsic part of the many-body wave 
function, one finds two different three-determinantal combinations, namely 

<M = I ITT) + e^^'/^l Tit) + e-'^'/^l TTi), (2.29) 

and 

Cr(7o) = I ITT) + e-^^'/^l Tit) + e'^'/^l rti), (2.30) 

where 70 = denotes the azimuthal angle of the vertex of the equilateral triangle 
associated with the original spin-down orbital in | itt)- We note that, unlike the 
intrinsic UHF Slater determinant, the intrinsic wave functions ^^t,. and ^^1^ here are 
eigenstates of the square of the total spin operator (S = with quantum 

number s = 1/2. This can be verified directly by applying to themj^ 

To restore the circular symmetry in the case of a (0,N) ring arrangement, one applies 



the projection operator (2.19). Note that the operator Vl is a direct generalization of 
the projection operators for finite point-groups discussed in Ref. |25] to the case of the 
continuous cyclic group Coo [the phases exp(i7L) are the characters of Coo]- 

The symmetry- restored projected wave function, \E'prj, (having both good total 
spin and angular momentum quantum numbers) is of the form, 

2n^PKj= /'''d7$L(7)e'^^ (2.31) 



where now the intrinsic wave function [given by (2.29) or (2.30)] has an arbitrary 
azimuthal orientation 7. We note that, unlike the phenomenological Eckardt-frame 
model |128t I132j where only a single product term is involved, the PRJ wave function in 



(2.31) is an average over all azimuthal directions over an infinite set of product terms. 
These terms are formed by multiplying the intrinsic part ^^tri'^) with the external 
rotational wave function exp(z7L) (the latter is properly characterized as "external", 
since it is an eigenfunction of the total angular momentum L and depends exclusively 
on the azimuthal coordinate 7)]^ 

The operator i?(27r/3) = exp(— i27rI//3) can be applied onto ^E'prj in two different 
ways, namely either on the intrinsic part or the external part exp(i7L). Using 
(2.29) and the property R{27r/3)^f^,^ = exp(-27ri/3)$j^j„ one finds, 

J?(27r/3)^PRj = exp(-27ri/3)^pRj, (2.32) 

% For the appropriate expression of S^, see equation (6) in Ref. |33]. 

Although the wave functions of the Eckardt-frame model are inaccurate compared to the PRJ ones 



[see (2.31)], they are able to yield the proper magic angular momenta for (0, A^) rings. This result is 
intuitively built in this model from the very beginning via the phenomenological assumption that the 
intrinsic wave function, which is never specified, exhibits Cjyy point-group symmetries. 
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Figure 7. (Color online) The UHF case exhibiting breaking of the circular symmetry 
for = 3 and Sz = 3/2 at Rw = 10 and B = 0. (a-c): real orbitals. (d): 
the corresponding electron density (ED). The choice of the remaining parameters is: 
hujQ — 5 meV and m* = 0.067me. Distances are in nanometers. The real orbitals are 
in 10^^ nm~^ and the total electron density in 10"'* nm~^. The arrows indicate the 
spin direction. 



from the first alternative, and 

J?(27r/3)*PRj = exp(-27rLi/3)^pRj, (2.33) 

from the second alternative. Now if \l/pRj 7^ 0, the only way that equations ( 2.32[ ) and 
(2.33) can be simultaneously true is if the condition exp[27r(L — l)i/3] = 1 is fulfilled. 
This leads to a first sequence of magic angular momenta associated with total spin 
s = 1/2, i.e., 

L = 3k + 1, k = 0, ±1, ±2, ±3, . . . (2.34) 



Using (2.30) for the intrinsic wave function, and following similar steps, one can 
derive a second sequence of magic angular momenta associated with good total spin 
s = 1/2, i.e., 

L = 3A; - 1, A; = 0, ±1, ±2, ±3, . . . (2.35) 

In the fully spin-polarized case, the UHF determinant is portrayed in figure [7] 
This UHF determinant, which we denote as | |tt ), is already an eigenstate of with 
quantum number s = 3/2. Thus only the rotational symmetry needs to be restored, that 
is, the intrinsic wave function is simply $J^tr(7o) = I TTT )• Since -R(27r/3)$j^tj. = $• 



intr' 



the condition for the allowed angular momenta is exp[— 27rLi/3] = 1, which yields the 
following magic angular momenta, 

L = 3k, k = 0, ±1, ±2, ±3, . . . (2.36) 
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We note that in high magnetic fields only the fully polarized case is relevant and 



that only angular momenta with k > enter in (2.36) (see Ref. [2l])- In this case, 
in the thermodynamic limit, the partial sequence with A; = 2g + 1, g = 0, 1, 2, 3, . . . is 
directly related to the odd filling factors u = l/(2g + 1) of the fractional quantum Hall 
effect [via the relation u = N{N — 1)/(2L)]. This suggests that the observed hierarchy 
of fractional filling factors in the quantum Hall effect may be viewed as a signature 
originating from the point group symmetries of the intrinsic wave function ^intr, and 
thus it is a manifestation of symmetry breaking at the UHF mean-field level. 

We further note that the discrete rotational (and more generally rovibrational) 
collective spectra associated with symmetry-breaking in a QD may be viewed as 
finite analogs to the Goldstone modes accompanying symmetry breaking transitions 
in extended media (see Ref. [ID])- Recently there has been some interest in studying 
Goldstone-mode analogs in the framework of symmetry breaking in trapped BECs with 
attractive interactions [88] . 

2. 3. The symmetry breaking dilemma and density functional theory 

Density functional theory (and its extension for cases with a magnetic field known 
as current density functional theory) was initially considered |3] (and was extensively 
applied [3l[92l|9l]) as a promising method for studying 2D semiconductor QDs. However, 
it soon became apparent [221 1231 |25l |16] that density functional approaches exhibited 
severe drawbacks when applied to the regime of strong correlations in QDs, where the 
underlying physics is associated with symmetry breaking leading to electron localization 
and formation of Wigner molecules. The inadequacies of density functional approaches 
in the field of QDs have by now gained broad recognition [HI I134[ 1135] . 

In particular, unlike the Hartree-Fock case for which a consistent theory for 
the restoration of broken symmetries has been developed (see, e.g., the earlier Refs. 
[18] |30| [3T] |32| [33] : for developments in the area of quantum dots, see the more recent 
Refs. [221 123, 1211 [23 US]), the breaking of space symmetry within the spin-dependent 
density functional theory poses [136] a serious dilemma. This dilemma has not been 
resolved [137] to date; several remedies are being proposed, but none of them appears to 
be completely devoid of inconsistencies. In particular, a theory for symmetry restoration 
of broken-symmetry solutions [134^ 1135] within the framework of density functional 
theory has not been developed as yet. This puts the density functional methods in a clear 
disadvantage with regard to the modern fields of quantum information and quantum 



computing; for example, the description of quantum entanglement (see section 5.1.4 
below) requires the ability to calculate many-body wave functions exhibiting good 
quantum numbers, and thus it lies beyond the reach of density functional theory. 

Moreover, due to the unphysical self-interaction error, the density-functional theory 
becomes erroneously more resistant to space symmetry breaking [138] compared to the 
UHF (which is free from such an error), and thus it fails to describe a whole class 
of broken symmetries involving electron localization, e.g., the formation at B = 
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of Wigner molecules in quantum dots [20l SB] and in thin quantum wires |139] . the 
hole trapping at Al impurities in silica |140] . or the interaction driven localization- 
delocalization transition in d- and /- electron systems, like Plutonium |141j . 

Recently, the shortcomings of the density functional theory to properly describe 
magnetic phenomena (such as exchange coupling constants associated with symmetry 
breaking of the total spin) has attracted significant attention in the quantum chemistry 
literature (see, e.g., Refs. [Ti^ [TiHl [Tii] ) . 

2.4- More on symmetry restoration methods 

In the framework of post-Hartree-Fock hierarchical approximations, projection 
techniques are one of methods used to treat correlations beyond the unrestricted 
Hartree-Fock. Two other methods are briefly discussed in this section, i.e., the method 
of symmetry restoration via random phase approximation (RPA) and the generator 
coordinate method (GCM). 

2.4- 1- Symmetry restoration via random phase approximation. This method introduces 
energy correlations by considering the effect of the zero-point motion of normal 
vibrations associated with the small amplitude motion of the time-dependent-Hartree- 
Fock mean field (which is equivalent to the RPA). In the case of space symmetry 
breaking, one of the RPA vibrational frequencies vanishes, and the corresponding motion 
is associated with the rotation of the system as a whole (rotational Goldstone mode), 
with a moment of inertia given by the so-called Thouless-Valatin expression |145j . 

The method has been used to calculate correlation energies of atomic nuclei 
[146[ 1147] and most recently to restore the broken symmetry in circular quantum dots 
[148] (mainly for the case of two electrons at zero magnetic field) . As discussed in Ref . 
[148] . restoration of the total spin cannot be treated within RPA. 

2.4-2. The generator coordinate method. The projection techniques by themselves do 
not take into account quantum correlation effects arising from the vibrations and other 
large-amplitude intrinsic collective distortions of the Wigner molecule. For the inclusion 
of the effects of such collective motions, a natural extension beyond projection techniques 
is the generator coordinate method (see Ch 10 in Ref. (18j). Unlike the RPA, the GCM 
can treat large-amplitude collective motion in combination with the retoration of the 
total spin. Indeed, it has been shown that the RPA harmonic vibrations are a limiting 
small-amplitude case of the large-amplitude collective motion described via the generator 
coordinate method [TS] . 

The GCM represents an additional step in the hierarchy of approximations 
described in section 11.51 and its use will result in a further reduction of the difference 
from the exact solutions. The GCM is complicated and computationally more expensive 
compared to projection techniques. Recent computational advances, however, have 
allowed rather extensive applications of the method in nuclear physics (see, e.g., Ref. 
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[38]). As yet, applications of the GCM to quantum dots or trapped atomic gases have 
not been reported. 

The GCM employs a very general form for the trial many-body wave functions 
expressed as a continuous superposition of determinants (or permanents for 

bosons), i.e., 

1$) = / d[a]f[ama]), (2.37) 



where [a] = (ai,a2, . . . , a^) is a set of collective parameters depending on the physics 
of the system under consideration. An example of such parameters are the azimuthal 



angles (71, 72, . . . , 7^) in the REM trial wave function (2.25). Of course the crucial 



difference between the REM wave function (2.25) and the general GCM function (2.37) 



is the fact that the weight coefficients f[a] in the former are known in advance (they 
coincide with the characters of the underlying symmetry group), while in the latter they 
are calculated numerically via the Hill- Wheeler-Griffin equations |149[ 1150] 



d[a']h{a,a')f[a'] = E J d[a']n{a,a')f[a'], (2.38) 
where E are the eigenenergies, and 

h{a,a) = {^[a]\n\^[a]), (2.39) 

n{a,a) = {^[a]\^[a]) (2.40) 



are the Hamiltonian and overlap kernels. The Hill- Wheeler-Griffin equation (2.38) is 
usually solved numerically by discretization; then one can describe it as a diagonalization 
of the many-body hamiltonian in a nonorthogonal basis formed with the determinants 
Ma]). 

An example of a potential case for the application of the GCM is an anisotropic 
quamtum dot {( < 1, with ( = u^/uy) In this case, one cannot use projection 
techniques to restore the total angular momentum, since the external confinement 
does not possess circular symmetry. Application of the GCM, however, will produce 
numerical values for the expansion coefficients /[7], and these values will reduce to 
exp[iX]g=i IqLq] for the circular case C = 1 [while the GCM wave function will reduce to 



the REM wave function (2.25)]. It is apparent that the GCM many-body wave function 
changes continuously with varying anisotropy (, although the symmetry properties of 
the confinement potential change in an abrupt way at the point C 1- 



3. Symmetry breaking and subsequent symmetry restoration for neutral 
and charged bosons in confined geometries: Theoretical framework 

3.1. Symmetry breaking for bosons, Gross-Pitaevskii wave functions, and permanents 

Mean-field symmetry breaking for bosonic systems is transparent in the context 
of two-component condensates, where each species is necessarily associated with a 
different space orbital |15H 1152] . For one species of bosons, symmetry breaking can 
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be considered through a generahzation of the UHF method of different orhitals for 
different spins known from the case of electrons in quantum chemistry and in quantum 
dots (section 2.1.1). Indeed, as shown in Refs. [601 ESI I153j . one can allow each 
bosonic particle to occupy a different orbital 0i(ri). Then the permanent \^n) = 
Perm[0i(ri), . . . , 07v(i"Ar)] serves as the many-body wave function of an unrestricted 
Bose-Hartree-Fock (UBHF) approximation. This wave function reduces to the Gross- 
Pitaevskii form with the restriction that all bosons occupy the same orbital 0o(i")) i-e-, 
= nili0o(i"j), and 0o(i") is determined self-consistently at the restricted Bose- 
Hartree-Fock (RBHF) level via the equation |154j 

[Ho{r,) + {N-1)J dv2<Pl{v2)v{v,, r2)M^2)]M^i) = ^o</>o(ri). (3.1) 

Here f(ri,r2) is the two-body repulsive interaction, which is taken to be a contact 
potential, vs = gS{ri — r2), for neutral bosons, or the Coulomb repulsion vc = 
e^Z^/(K|ri — r2|) for charged bosons. The single-particle Hamiltonian is given by 

Ho{v) = -^^VV(2m) + (3.2) 

where ujq characterizes the circular harmonic confinement, and where we have considered 
a non-rotating trap. 

Compared to the fermionic case (see section 2.1.1 ), the self-consistent determination 
of the UBHF orbitals is rather complicated and numerically highly demanding (see 
section 2.5.3 in Ref. |155] ). In fact, the self-consistent UHBF equations cannot be 
put into a standard (canonical) eigenvalue-problem form due to two reasons: (i) unitary 
transformations cannot be used to simplify the equations, since the permanent of the 
product of two matrices does not factorize into a product of two simpler terms (unlike 
the electronic case where the determinant of the product of two matrices is equal to 
the product of the corresponding determinants) and (ii) as a result of boson statistics, 
the bosonic orbitals cannot be assumed to be (and remain) orthogonal, which leads to 
additional coupling terms between the non-orthogonal orbitals. 

In the literature |156] . an attempt has been made to derive unrestricted self- 
consistent equations for bosons by disregarding point (ii) mentioned above and invoking 
the assumption of orthonormal orbitals. Such equations of course are not of general 
validity, although they appear to be useful for describing fragmentation of Bose 
condensates in double wells. 

The difficulties of the self-consistent treatment can be bypassed and the UBHF 
problem can be simplified through consideration of explicit analytic expressions for the 
space orbitals 0i(rj). In particular, for repulsive interactions, the bosons must avoid 
each other in order to minimize their mutual repulsion, and thus, in analogy with the 
case of electrons in QDs, one can take all the orbitals to be of the form of displaced 
Gaussians, namely, 

0,(r,) = TT-i/V-i exp[-(r, - a,)V(2a2)]. (3.3) 

The positions aj describe the vertices of concentric regular polygons, with both the 
width a and the radius a = | | of the regular polygons determined variationally through 
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minimization of the total energy 

^UBHF = ($iv|H|$iv)/($w|$7v), (3.4) 

where 

N N 

n = J2Ho{r,) + J2vir,,v^) (3.5) 

i=l i<j 

is the many-body Hamiltonian. 

With the above choice of locahzed orbitals the unrestricted permanent \^n) breaks 
the continuous rotational symmetry. However, the resulting energy gain becomes 
substantial for stronger repulsion. Controlling this energy gain (the strength of 
correlations) is the ratio Rs (or Rw) between the strength of the repulsive potential 
and the zero-point kinetic energy. Specifically, for a 2D trap, one has 

Rs = gm/{2Txn^) (3.6) 



for a contact potential (for Rw-, see section 2.1.2). Note that in this section, we refer to 



the case of a non-rotating trap; the generalization to rotating traps is presented later in 
section |H1 

3.2. Repelling bosons in harmonic traps: Restoration of broken symmetry. 

Although the optimized UBHF permanent \^n) performs very well regarding the total 
energies of the trapped bosons, in particular in comparison to the resctricted wave 
functions (e.g., the GP anzatz), it is still incomplete. Indeed, due to its localized 
orbitals, \^n) does not preserve the circular (rotational) symmetry of the 2D many-body 
Hamiltonian H. Instead, it exhibits a lower point-group symmetry, i.e., a C2 symmetry 
for N = 2 and a C5 one for the (1,5) structure of = 6 (see section [s] below). As a 
result, |$Ar) does not have a good total angular momentum. In analogy with the case 
of electrons in quantum dots, this paradox is resolved through a post-Hartree-Fock step 
of restoration of broken symmetries via projection techniques [231 EH EHl [60], yielding 
a new wave function I^P^^) with a definite angular momentum L, that is 

2vr|^^^i)= ^d7|$;v(7))e'^^ (3.7) 
Jo 

where is the original UBHF permanent having each localized orbital rotated 

by an azimuthal angle 7, with L being the total angular momentum. The projection 
yields wave functions for a whole rotational band. Note that the projected wave function 



|\E'^^) in (3.7) may be regarded as a superposition of the rotated permanents |$jv(7)), 
thus corresponding to a "continuous-configuration-interaction" solution. 
The energies of the projected states are given by 

El^^ = (vl/^^^|7^|vl/^5)/(vI/^^^|v^^^^). (3.8) 
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4. Other many-body methods 

4.I. Exact diagonalization methods: Theoretical framework 

We will discuss the essential elements of the exact-diagonalization method here by 
considering the special, but most important case of a many-body Hilbert space defined 
via the restriction that the single-particle states belong exclusively to the lowest Landau 
level. For 2D electrons, this LLL restriction is appropriate in the case of very high B. 
For rotating bosons in a harmonic trap, the LLL restriction is appropriate for O ~ cuq 
and a very weak repulsive contact potential. The particulars of the EXD method for 
quantum dots in the case of field-free (and/or low B) conditions will be discussed below 
in section I5.1.2[ 

For sufficiently high magnetic field values (i.e., in the fractional quantum Hall effect, 
regime), the electrons are fully spin-polarized and the Zeeman term (not shown here) 
does not need to be considered. In the -B — > oo limit, the external confinement V{x, y) 
can be neglected, and the many-body H can be restricted to operate in the lowest 
Landau level, reducing to the form [211 [261 IS2l [53] 

HLLL = iV^ + EE;|r> (4.1) 

where Uc = eB/{m*c) is the cyclotron frequency. Namely, one needs to diagonalize the 
interaction Hamiltonian only. 



For the case of rotating bosons in the LLL, one needs to replace in (4.1 ) the Coulomb 
interaction by gS{ri — rj) and the cyclotron frequency by 2Q (see the Appendix for the 
details of the equivalence between magnetic field B and rotational frequency Q; in short 
2Q). 

For a given total angular momentum L = Y^^=i ^k, the exact-diagonalization A^- 
body wave function is a linear superposition of Slater determinants for fermions (or 
permanents for bosons) "^{J) made out of lowest-Landau-level single-particle wave 
functions (see the Appendix), 

Mz) = -^(^)'e-'^^'^^\ (4-2) 



AVvf/I VA 



where A = J2hc/{eB) = IbV^ for the case of electrons in QDs {Ib being the magnetic 



length), and A = yh/{muJo) for the case of bosons in rotating traps. In (4.2), we used 
complex coordinates z = x + iy, instead of the usual vector positions r = (x, y); below 
we will use either notation interchangeably as needed. 

Thus, the many-body EXD wave function is written as 

$^^°(;^l,;^2,...,^^) = ECjV[/(J) (4.3) 

J=l 

with the index J denoting any set of A^ single-particle angular momenta {/i, I2, . . . ,In} 
such that 

h < k < . . . < In (4.4) 



Symmetry breaking and quantum correlations 



37 



for fermions and 

h<l2<...<lN (4.5) 

for bosons, the absence or presence of the equal signs being determined by the different 
statistics between fermions and bosons, respectively. 

Using the expansion (4.3), one transforms the many-body Schrodinger equation 

Z2,..., z^) = z,,..., z^), (4.6) 

into a matrix eigenvalue problem. When the single-particle states are orthonormal 
[like the LLL ones in (4.2)], the matrix elements {\E'(/)|7i|\E'(J)) between two Slater 
determinants are calculated using the so-called Slater rules (see, e.g., Ch 2.3.3 of 
Ref. [in]). For the case of bosons, the correpsonding rules for the matrix elements 
(\E'(/)|?i|^(J)) between two permanents are given in the appendix of Ref. [72] . 

We remark here that the calculation of energies associated with projected wave 
functions [see, e.g.. Equation (3.8)] requires calculation of similar matrix elements 
between two Slater determinants (or permanents) with non-orthogonal orbitals; the 
corresponding formulas for the case of fermions can be found in Ch 6.3. of Ref. |157] . 
and for the case of bosons in Ref. [155] . 

Of course a necessary ingredient for the application of the above rules is the 
knowledge of the matrix elements of the two-body interaction v{ri, r2), i.e., 

Vap,A = J dridr2C(l)0;(2)t;(l,2)0^(l)0A(2). (4.7) 

In the general case, these two-body matrix elements need to be calculated 
numerically. For the simpler case specified by the LLL orbitals (4.2), the two-body 
matrix elements are given by analytic expressions. In particular, for the Coulomb 
repulsion, see Refs. ^2E[ I158j : for a contact potential, see Ref. [159j . 

The dimension T) of the Hilbert space increases very fast with the number of 
particles N and the value of the total angular momentum L, and is is controlled by 
the maximum allowed single-particle angular momentum /max, such that Ik < /max, 
1 < A; < A^. By varying /max, we can check that this choice produces well converged 
numerical results. 

For the solution of the large scale, but sparse, EXD matrix eigenvalue problem 
associated with the special Hamiltonian ./^lll [or the general one in (2.1)], we use the 
ARPACK computer code pOj . 

The availabilty of analytic expressions for the two-body interaction has greatly 
facilitated exact-diagonalization calculations in the lowest Landau level (appropriate 
for quantum dots at high B), and in this case (starting with Refs. [261 1161] ) 
diagonalization of large matrices of dimensions of order 500,000 x 500,000 has become 
a commom occurrence. For circular quantum dots, similar analytic expressions for 
the matrix elements of the Coulomb interaction between general Darwin-Fock orbitals 
|162[ 1163] (i.e., the single-particle orbitals of a circular 2D harmonic oscillator under 
a perpendicular magnetic field B) are also available [1291 1164] . but they are not 
numerically as stable as Tsiper's expressions [158] in the lowest Landau level. 
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Exact-diagonalization calculations for field- free (and/or low B) conditions have 
been presented in several papers. Among them, we note the exact-diagonalization 
calculations of Refs. [771 EH MEi [M [Ml MB HBHl UTO]. EXD calculations employing 
Coulombic two-body matrix elements that are calculated numerically have also been 
reported for elliptic quantum dots (see section 5.1.2 below). Furthermore, some authors 
have used the method of hyperspherical harmonics [171] for circular quantum dots, 
while others have carried out exact-diagonalization calculations for quantum dots with 
a polygonal external confinement [172j . 

Concerning EXD calculations in the lowest Landau level, we mention Refs. 
[23 [SSI [ffl [M [M [Ml [nSl [na inSl HZS] for the case of quantum dots (high 
B) and Refs. [671 (Ml ISHl [721 1177] 1178] for the case of bosons in rapidly rotating traps. A 
version of EXD in the LLL uses a correlated basis constructed out of composite-fermion 
wave functions [176j . while another exact-diagonalization version used non-orthogonal 
floating Gaussians in the place of the usual single-particle states (4.2) in the LLL. 

For two electrons in a single quantum dot, exact calculations have been carried 
out through separation into center-of-mass and relative coordinates [50l 11791 1180] . In 
addition, EXD calculations have been reported for two electrons in a double quantum 
dot [891IT8T]. 

It is of interest to note that the EXD approach is also used in other fields, but 
under different names. In particular, the term "shell model calculations" is used in 
nuclear theory, while the term "full configuration interaction" is employed in quantum 
chemistry. 



4.I.I. An example involving spin-resolved CPDs. Here we present an example of an 
EXD calculation exhibiting formation of a Wigner molecule in quantum dots. The case 
we chose is that of = 3 electrons under zero magnetic field in an anisotropic quantum 
dot with huJx = 4.23 meV and huy = 5.84 meV (i.e., with an intermediate anisotropy 
( = uj^/ujy = 0.724) and dielectric constant k = 1 (strong interelectron repulsion). In 
particular, figure [8] displays results for the ground state of the three electrons with total 
spin s = 1/2 and spin projection 5*2 = 1/2. 

The electron density in figure [8|a) has the shape of a diamond and suggests 
formation of a Wigner molecule resonating between two isosceles triangular isomers 
(which are the mirror image of each other). The detailed interlocking of the two 
triangular configurations is further revealed in the spin-resolved CPD that is displayed in 
figure [Sj^b). It can be concluded that one triangle is formed by the points Ri = (0, —20) 
nm, R2 = (—43, 10) nm, and R3 = (43, 10) nm, while the second one (its mirror with 
respect to the x-axis) is formed by the points R'l = (0,20) nm, R'2 = (—43, —10) nm, 
and R'3 = (43, -10) nm. 

The two-triangle configuration discussed for three electrons above may be seen 
as the embryonic precursor of a quasilinear structure of two intertwined "zig-zag" 
crystalline chains. Such double zig-zag crystaline chains may also be related to the 
single zig-zag Wigner-crystal chains discussed recently in relation to spontaneous spin 
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Figure 8. (Color online) Exact-diagonalization electron densities (EDs) and spin- 
resolved CPDs for N = 3 electrons in an anisotropic quantum dot at zero magnetic field 
{B = 0) and in the case of a strong Coulomb repulsion with a dielectric constant k = I. 
(a-b) Results for the ground state (which has total spin s = 1/2 and spin projection 
Sz = 1/2) when Tiojx = 4.23 meV and Tiujy = 5.84 meV (i.e., for an intermediate 
anisotropy ^ = 0.724). (c-f) Results for the first excited state (which has also total 
spin s = 1/2 and spin projection Sz — 1/2) when hujx — 3.137 meV and hujy — 6.274, 
meV (i.e., for a strong anisotropy C = 0.5). The thick arrows and solid dots in the 



CPDs indicate the spin direction ctq and position Tq of the fixed electron [see ( 1.2 )]. The 
thin arrows indicate the spin direction of the remaining two electrons. The effective 
mass is m* = 0.070TOe for the intermediate anisotropy (a-b) and m* = 0.067me for the 
strong anisotropy (c-f). Lengths are in nanometers. The vertical axes are in arbitrary 
units. 



polarization in quantum wires |182l I183j . 

For strong anisotropics (e.g., ( < 1/2), the three electrons form a straightforward 
linear Wigner molecule [see the electron density in figure [8|^c)], and the spin- resolved 
CPDs can be used to demonstrate |184j formation of prototypical entangled states, like 
the so-called W states |185l 1186] . From the CPDs [displayed in figure [sj^d), figure [sj^e) 
and figure |8](f)] of the first excited state (having s = 1/2 and 5*^ = 1/2), one can infer 
that its intrinsic spin structure is of the form |tTi)~|iTT)- The ground state (not 
shown) of this linear Wigner molecule has also a total spin s = 1/2 and spin projection 
Sz = 1/2, and its intrinsic spin structure corresponds to a form 2| llf ) — | fH ) — | ||| ) 

4-2. Particle localization in Monte- Carlo approaches 

Quantum Monte-Carlo (MC) approaches |187] have been successfully used in many 
areas of condensed-matter physics; they are divided in two main branches, path-integral 
MC (PIMC) and variational/diffusion MC (V/DMC). Unlike the exact diagonahzation, 
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quantum MC approaches cannot calculate excited states, and they are restricted to the 
description of ground-state properties 

Applied to the case of circular quantum dots at zero magnetic field, PIMC 
calculationss |188[ 11891 11901 119H 1192] have been able to reproduce and describe electron 
localization with increasing and formation of Wigner molecules. In 2D quantum 
dots, a focus of the PIMC studies |188t 1191] has been the determination of the critical 
value, R^, for the Wigner parameter at which the phase transition from a Fermi liquid 
to a Wigner molecule occurs. Naturally, only an estimate of this critical value can be 
determined, since the phase transition is not sharp, but smooth, due to the finite size of 
the quantum dot. Obtaining a precise value of is also hampered by the variety of 
criteria employed by different researchers in the determination of this transition (e.g., 
height of localized density humps, appearance of a hump at the center of the dot, etc.). 

In the literature of PIMC studies |188t 1191] , one finds the critical value R^r ~ 4, 
which is in agreement with exact-diagonalization studies [77]. This is also in general 
agreement with the estimate i?^ ~ 1 based on the abrupt onset of spatial symmetry 
breaking in unrestricted Hartree-Fock [20j. Of course the unrestricted-Hartree-Fock 
estimate has to be refined through the subsequent step of symmetry restoration. We 
believe that it is most appropriate to consider these two estimates mentioned above 
as the lower and upper limit of a transition region. The important conclusion is that 
the transition to Wigner crystallization in quantum dots takes place for much higher 
electron densities compared to the infinite two-dimensional electron gas (for which a 
value R^ ~ 37 [193] has been reported.)]^ 

A disadvantage of the PIMC method is that the case of an applied magnetic field 
cannot be easily incorporated in its formalism, and therefore related studies have not 
been reported. Other well known difficulties are the fermion sign problem and the 
nonconservation of total spin [TTJ 1188] . 

Commenting on the other main branch of quantum Monte Carlo, i.e., the 
variational/diffusion MC, we wish to stress the crucial role played by the general 
form of the trial wave function used. Indeed, an early V/DMC study [194] using a 
single configurational state function (i.e., a primitive combination of products of Slater 
determinants for the two spin directions that is an eigenstate of the total angular 
momentum L, the square of the total spin S^, and the total-spin projection Sz) was 
unable to describe the formation of Wigner molecules in quatum dots at zero magnetic 
field. Another V/DMC study [195] managed to demonstrate electron localization, but 
at the cost of using a single product of two Slater determinants (multiplied by a Jastrow 
factor) which violated the conservation of both the total angular momentum and total 
spin (without the possibility of further corrections related to symmetry restoration). 

Most recently, more sophisticated trial wave functions involving a large number 
of configurational state functions with good total angular momentum and total spin 

* Often the Wigner-Seitz radius r^, in units of the effective Bohr radius a*g ~ fi"^ k / {m* e^) of the 
quantum dot, is used instead of the Wigner parameter Rw (denoted some times by A). In these units, 
one has Ts ~ Rw- 
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have been employed, which enabled eventually confirmation of the formation of Wigner 
molecules via V/DMC methods, both at zero |78j and high magnetic field |196j . 

There are, however, disagreements between the V/DMC results |197j and those 
from PIMC and EXD calculations concerning the details of Wigner-molecule formation 
in circular quantum dots in the absence of an applied magnetic field. In particular, 
these disagreements focus on the density scale for the cross-over and the strength of 
azimuthal and radial electron correlations as a function of Rw- 

Such disagreements remain an open question for two reasons: 

(i) The criterion of lowest energy (evoked by the V/DMC approaches) is not sufficient 
to guarantee the quality of the variational many-body wave function. A counterexample 
to this lowest-energy criterion was presented by us for the case of the Laughlin wave 
functions in Refs. [2H E2] (see also section 6.4). Most recently, this point was also 
illustrated within the framework of variational Monte Carlo calculations |198] . 

(ii) The V/DMC studies for larger N [THl I197j have presented only calculations for 
CPDs. However, due to the presence of dummy integrations in (1.1 ) (which result in an 
averaging over the remaining N — 2 particles), the ability of the CPDs to portray the 
intrinsic crystalline structure of the Wigner molecule diminishes with increasing A^. As 
a result, higher-order correlation functions, like A^-point correlations, may be required. 
The fact that higher-order correlation functions reflect the crystalline correlations more 
accurately than the CPDs was illustrated for the case of rotating boson molecules in 



Ref. [72] (see also section 8.2). 

A detailed comparison between ground-state energies calculated with quantum MC 
and exact-diagonalization methods can be found in Ref. [77j. For a comparison between 
variation-before-projection (see section 2.2) and V/DMC total energies, see Ref. |122] . 



5. The strongly correlated regime in two-dimensional quantum dots: The 
two-electron problem and its significance 

In section 2 and section 3, we focused on the general principles and the essential 
theoretical framework of the method of symmetry breaking and of subsequent symmetry 
restoration for finite condensed-matter systems. In addition, in section 4, we presented 
the basic elements of the exact-diagonalization approach. In the following four sections, 
we will focus on specific applications and predictions from these methods in the 
field of semiconductor quantum dots and of ultracold bosons in harmonic traps, in 
particular regarding the emergence and properties of Wigner molecules under various 
circumstances. At the same time we will continue to elaborate and further expand on 
more technical aspects of these methods. 

In this section, we start by concentrating on the description of two-electron 
molecules in QDs. A discussion on the importance of the two-electron problem is given 
in section 15.31 
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5.1. Two-electron elliptic dot at low magnetic fields 

Here, we present an exact diagonalization and an approximate (generalized Heitler- 
London, GHL) microscopic treatment for two electrons in a single elliptic QD specified 
by parameters that correspond to a recently fabricated experimental device [H]. 
The two-dimensional Hamiltonian for the two interacting electrons is given by 

n = H{r^) + Hir^) + le^it^r.^), (5.1) 

where the last term is the Coulomb repulsion, k (12.5 for GaAs) is the dielectric constant, 
and ri2 = |ri — r2|. The prefactor 7 accounts for the reduction of the Coulomb strength 
due to the finite thickness of the electron layer in the z direction and for any additional 
screening effects due to the gate electrons. H{r) is the single-particle Hamiltonian for an 
electron in an external perpendicular magnetic field B and an appropriate confinement 



potential (2.2). For an elliptic QD, the external potential is written as 

V{x,y) = ^m*{ujy + ujy). (5.2) 

Here the effective mass is taken to be m* = 0.07mo. In the Hamiltonian (2.2), we 
neglect the Zeeman contribution due to the negligible value {g* ~ 0) of the effective 
Lande factor in our sample [199] . 

5.1.1. Generalized Heitler- London approach. The GHL method for solving the 



Hamiltoninian (5.1) consists of two steps. In the first step, we solve selfconsistently 
the ensuing unrestricted Hartree-Fock equations allowing for lifting of the double- 
occupancy requirement (imposing this requirement gives the restricted HF method, 
RHF). For the 3^ = solution, this step produces two single-electron orbitals uifiij) 
that are localized left (L) and right {R) of the center of the QD [unlike the RHF 
method that gives a single doubly-occupied elliptic (and symmetric about the origin) 
orbital]. At this step, the many-body wave function is a single Slater determinant 
^uhf(1 i) = |ml(1 ])ur{2 I)) made out of the two occupied UHF spin-orbitals 
■Ui(l 1) = ■Ui(ri)a(l) and ur{2 [) = uji{r2)P(2), where a(/3) denotes the up (down) 
[T (i)] spin. This UHF determinant is an eigenfunction of the projection 5*^ of the total 
spin S = Si + §2, but not of 5*^ (or the parity space- reflection operator). 

In the second step, we restore the broken parity and total-spin symmetries by 



applying to the UHF determinant the projection operator (2.18). For two electrons 



this operator reduces to Pgpin = 1 =F ^12, where the operator ti7i2 interchanges the spins 



spin 



of the two electrons; the upper (minus) sign corresponds to the singlet. The final result 
is a generalized Heitler-London two-electron wave function \E'Q^L(ri, r2) for the ground- 
state singlet (index s) and first-excited triplet (index t), which uses the UHF localized 
orbitals, 

^GHL(ri, r2) oc (ML(ri)Mi?(r2) ± ML(r2)MR(ri))x''*, (5.3) 

where x'^'* = ('^(l)/5(2) =F ck^(2)/9(l)) is the spin function for the 2e singlet and triplet 
states. The general formalism of the 2D UHF equations and of the subsequent 
restoration of broken spin symmetries was presented in section 2.2. 
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The use of optimized UHF orbitals in the generahzed Heitler-London method is 
suitable for treating single elongated QDs |16], including the special case of elliptically 
deformed ones discussed in this section. The GHL is equally applicable to double QDs 
with arbitrary interdot-tunneling coupling [36]. In contrast, the Heitler-London (HL) 
treatment [S2] (known also as the simple Valence bond) , where non-optimized "atomic" 
orbitals of two isolated QDs are used, is appropriate only for the case of a double dot 
with small interdot-tunneling coupling 

The orbitals ul,r{tc) are expanded in a real Cartesian harmonic-oscillator basis, i.e., 

«L,.j(r) = f:cf'V(r), (5.4) 
i=i 

where the index j = {m,n) and v^j(r) = Xm{x)Yn{y), with Xm{Yn) being the 
eigenfunctions of the one-dimensional oscillator in the x{y) direction with frequency 
Ux{ujy). The parity operator V yields VXm{x) = {—l)"^Xm{x), and similarly for Yn{y). 
The expansion coefficients C^'^ are real for B = and complex for finite B. In the 
calculations we use K = 54 and/or K = 79, yielding convergent results. 



5.1.2. Exact diagonalization. In the EXD method, the many-body wave function 
is written as a linear superposition over the basis of non-interacting two-electron 
determinants, i.e., 

2K 

^^XD(ri,r2) =E^^V(l;0^(2;j)), (5.5) 

i<j 

where V(1;0 = <^i(l T) if 1 < « < and V(l; = <^i-A'(l I) if K + 1 < i < 2K 
[and similarly for ^jj{2,j)]. The total energies -Eexd coefficients are 

obtained through a "brute force" diagonalization of the matrix eigenvalue equation 
corresponding to the Hamiltonian in (5.1). The exact-diagonalization wave function 
does not immediately reveal any particular form, although, our calculations below show 
that it can be approximated by a GHL wave function in the case of the elliptic dot 
under consideration. 



5.1.3. Results and comparison with measurements. To model the experimental 
quantum dot device, we take, following Ref. [41], hojx = 4.23 meV, TiUy = 5.84 meV, 
K = 12.5, and 7 = 0.862. The corresponding anisotropy is Uy/ux = 1.38, indicating that 
the quantum dot considered here is closer to being circular than in other experimental 
systems [33 EO] . 

As shown in Ref. [31], the experimental findings can be quantitatively interpreted 
by comparing to the results of the EXD calculations for two electrons in an anisotropic 
harmonic confinement potential with the parameters listed above. All the states 
observed in the measured spectra (as a function of the magnetic field) can be 
unambiguously identified [31] with calculated ground-state and excited states of the 



two-electron Hamiltonian (compare figure [9] and figure 10). 
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magnetic field (T) 

Figure 9. (Color online) Differentiated current dl/dVpg at Vbias = 2.5 mV (the 
subscript pg denotes the plunger gate). Gray striped regions (red online) marked 
by symbols correspond to positive (peaks) dl/dVpg. The dark black region (also 
black online) corresponds to negative dl/dVpg. Electron numbers N are indicated. 
Transitions between the one-electron ground state and the 2e spin-singlet ground 
(excited for B > 3.8 T) state (Sq), spin-triplet excited (ground for B > 3.8 T) state 
(T+), spin-singlet excited state (S'2), and spin-triplet plus center-of-mass excited state 
(T+.cm) are labeled. 




B(T) 



Figure 10. (Color online) Calculated exact-diagonalization energy spectrum in a 
magnetic field, referenced to 2hy/Lol + ojI/A, of a 2e dot with anisotropic harmonic 
confinement (for the dot parameters, see text). We have adopted the notation 
{Nx,Ny,n,m), where {Nx,Ny) refer to the CM motion along the x and y axes 
and (n, m) refer to the number of radial nodes and angular momentum of the 
relative motion in the corresponding circular dot. Inset: The EXD spectrum of the 
corresponding circular dot. Only the (71, m) indices are shown, since iV^, — Ny = for 
all the plotted curves. Solid hues denote singlets. Dashed hues denote triplets. 
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Figure 11. (Color online) Comparison between the lowest-triplet/lowest-singlet 
energy splitting [J{B)] calculated with different methods and the experimental results 
(open squares). Solid line (online magenta): EXD. Dotted line (online green): GHL. 
Dashed line (online red): RHF. For the parameters used in the calculation to model 
the anisotropic QD, see text. 



Moreover, the calculated magnetic- field-dependent energy splitting, Jexd(-B) = 
-^exd(-^) ~ -^exd(-^)' between the two lowest singlet (Sq) and triplet (T+) states is 
found to be in remarkable agreement with the experiment [see figure [iT] . 

A deeper understanding of the structure of the many-body wave function can 
be acquired by comparing the measured J{B) with that calculated within the GHL 
and RHF approximations. To facilitate the comparisons, the calculated Jghl(-B) and 
</rhf(-B) curves are plotted also in figure 11, along with the exact-diagonalization result 
and the measurements. Both the RHF and GHL schemes are appealing intuitively, 
because they minimize the total energy using single-particle orbitals. It is evident, 
however, from figure 11 that the RHF method, which assumes that both electrons occupy 
a common single-particle orbital, is not able to reproduce the experimental findings. On 
the contrary, the generalized Heitler-London approach, which allows the two electrons to 
occupy two spatially separated orbitals, appears to be a good approximation. Plotting 
the two GHL orbitals [see figure 12 for the singlet state clearly demonstrates that the 
two electrons do not occupy the same orbital, but rather fill states that are significantly 
spatially separated. 

The UHF orbitals from which the GHL singlet state is constructed [see (5.3)] are 
displayed on figure 12 for both the B = and B = 3.8 T cases. The spatial shrinking of 
these orbitals at the higher 5-value illustrates the "dissociation" of the electron dimer 
with increasing magnetic field. The asymptotic convergence (beyond the ST point) of 
the energies of the singlet and triplet states, [i.e., J{B) — *• as S — *• oo] is a refiection of 
the dissociation of the 2e molecule, since the ground-state energy of two fully spatially 
separated electrons (zero overlap) does not depend on the total spin. We stress again 
that the RHF, which corresponds to the more familiar physical picture of a QD-Helium 
atom, fails to describe this dissociation, because J^up^B) diverges as the value of the 
magnetic field increases. 

In contrast to the RHF, the GHL wave function is able to capture the importance of 
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Figure 12. (Color online) Single-particle UHF orbitals (modulus square) that are 
used in the construction of the GHL wave function in ( |5.3| . Arrows indicate up and 
down spins. For the parameters used in the calculation to model the anisotropic QD, 
see text. Lengths in nm and orbital densities in 10~^ nm^^. 

B=0 CPDs B=3.8 T 




Figure 13. (Color online) CPDs extracted from the exact-diagonalization wave 
function for the singlet state for B — and B — 3.8 T. The CPD expresses the 
conditional probability for finding the second electron at position r given that the first 
electron is located at ro (denoted by a heavy solid dot). For the parameters used in 
the calculation to model the anisotropic QD, see text. Lengths in nm and CPDs in 
arbitrary units. 



correlation effects. Further insight into the inportance of correlations in our QD device 
can be gained through inspection of the conditional probability distributions [see 
(1.1 )] associated with the EXD solutions; see an illustration in figure 13 Indeed, already 
at zero magnetic field, the calculated CPDs provide further support of the physical 
picture of two localized electrons forming a state resembling an H2-type 
Wigner molecule [201 1188] . 



5.1.4- Degree of entanglement. Further connections between the strong correlations 
found in our microscopic treatment and the theory of quantum computing [48J can be 
made through specification of the degree of entanglement between the two localized 
electrons in the molecular dimer. For two electrons, we can quantify the degree of 
entanglement by calculating a well-known measure of entanglement such as the von 
Neumann entropy |l2l 1200] for indistinguishable particles. To this effect, one needs 
to bring the EXD wave function into a diagonal form (the socalled "canonical form" 
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Figure 14. Von Neumann entropy for the lowest singlet EXD state of the elliptic dot 
as a function of the magnetic field B. On the top, we show histograms for the jz^p 
coefficients [see of the singlet state at B = 3.8 T (left) and B = 8.0 T (right) 

illustrating the dominance of two determinantal configurations (in agreement with the 
generalized Heitler-London picture). Note the small third coefficient jzap = 0.081 for 
B — 8.0 T. For the parameters used to model the experimental device, see text. 



[2001 [201]), i.e., 

M 

nini^u r2) = E 2A: - 1)$(2; 2k)), (5.6) 

k=l 

with the $(i)'s being appropriate spin orbitals resulting from a unitary transformation 
of the basis spin orbitals ip{j)^s [see (5.5)]; only terms with 2;^ 7^ contribute. The 
upper bound M can be smaller (but not larger) than K (the dimension of the single- 
particle basis); M is referred to as the Slater rank. One obtains the coefficients of 
the canonical expansion from the fact that the Iz^l"^ are eigenvalues of the hermitian 
matrix A'^A [A, see (5.5), is antisymmetric]. The von Neumann entropy is given by 
S = — Y^kLi kfcP log2(|-2fcn with the normalization Y^kLi kfcP = 1- 

The EXD singlet has obviously a Slater rank M > 2. The von Neumann entropy for 
the EXD singlet (iS^xd) is displayed in figure 14 It is remarkable that iS|;xd increases 
with increasing B, but remains close to unity for large B, although the maximum allowed 
mathematical value is log2(-ft') [for example, for K = 79, log2(79) = 6.3]. The saturation 
of the entropy for large i? to a value close to unity refiects the dominant (and roughly 
equal at large B) weight of two configurations in the canonical expansion [see (5.6)] of 
the exact-diagonalization wave function, which are related [12] to the two terms in the 
canonical expansion of the GHL singlet. This is illustrated by the histograms of the 



Izll"^ coefficients for B = 3.8 T and B = 8.0 T in figure 14 (top). Notice that the ratio 
k2p/kiP refiects the extent of the overlap between the two GHL orbitals |12], with the 
ratio increasing for smaller overlaps (corresponding to a more complete dissociation of 
the Wigner molecule). 

The above discussion illustrates that microscopic calculations that are shown to 
reproduce experimental spectra [H] can be used to extract valuable information that 
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Figure 15. (Color online) The calculated spectrum of a two-electron circular 
parabolic quantum dot, with Rw — 200. The quantum numbers are (N, M, n, m) 
with N corresponding to the number of radial nodes in the center of mass (CM) 
wavefunction, and M is the CM azimuthal quantum number. The integers n and 
m are the corresponding quantum numbers for the electrons' relative motion (RM) 
and the total energy is given by EMM,nm — -E'^m + E^'^{n,\'m\). The spectrum 
may be summarized by the "spectral rule" given in the figure, with an efl^ective rigid 
moment of inertia C = 0.037 (corresponding to an angular momentum C = Tim), 
the phonon for the stretching vibration ttiOa — 3.50, and the phonon for the bending 
vibration coincides with that of the CM motion, i.e., fujJb — Tiuq — 2. The quantum 
numbers (Nq, Mo,no,m) specifying each rotational band are given at the bottom, 
with m = 0, 1, 2, . . . (the levels m — and to = 1 in each band may not be resolved 
on the scale of the figure). We note that the energy separation between levels in a 
given rotational band increases as (2to + 1) with increasing to, which is a behavior 
characteristic of a rigid rotor. All energies are in units of fiuJo/2, where ojq is the 
parabolic confinement frequency. 



allows assessment of the suitability of a given device for quantum computations. 

5.2. Two-electron circular dots at zero magnetic field 



In section |2.2[ we illustrated the formation of "rotating electron molecules" in the case 
of a two-electron circular QD, where one needs to consider restoration of the rotational 
symmetry as well, in addition to the restoration of the total spin. There, we focused on 
properties of the ground state (L = 0). 

In this section, we will further examine the excitation spectra of a two-electron QD 
by using the rather simple exact solution of this problem provided through separation 
of the center-of-mass and inter-electron relative-distance degrees of freedom [50]. The 



spectrum obtained for Rw = 200 (figure 15), exhibits features that are characteristic 



of a collective rovibrational dynamics, akin to that of a natural "near-rigid" triatomic 
linear molecule with an infinitely heavy middle particle representing the center of mass 
of the dot. This spectrum transforms to that of a "floppy" molecule for smaller value 
of Rw (i.e., for stronger confinements characterized by a larger value of ujq, and/or 
for weaker inter-electron repulsion), ultimately converging to the independent-particle 
picture associated with the circular central mean- field of the QD. 

Further evidence for the formation of the electron molecule and the emergence 
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Figure 16. CPD of the excited multi-vibrational state (1,0,1,0) of a 2e circular 
parabolic QD with Rw — 200 (see text for details) . The solid dot portrays the position 
of the reference point Tq = {do,0), where do = 2.6 is half the interelectron distance at 
the ground state (0,0,0,0). Distances are in units of loV2; the scale of the vertical 
axis is arbitrary. 



of a rovibrational spectrum was found through examination ^50j of the conditional 
probabihty distributions for various states (A^, M, n, m) (see the caption of figure 15 for 
the precise meaning of these quantum numbers labehng the spectra). As an example, 
we display in figure 16 the CPD for the bottom state (m = 0) of the rotational band 
(1, 0, 1, m) (not shown in figure 15 ); it reveals that this state corresponds to a vibrational 
motion of the electron molecule both along the interelectron axis (one excited stretching- 
mode phonon; see figure 15 ) and perpendicularly to this axis (two excited bending-mode 
phonons; see figure 15). 

It is instructive to note here certain similarities between the formation of a 
"two-electron molecule" in man-made quantum dots, and the collective (rovibrational) 
features observed in the electronic spectrum of doubly-excited helium atoms |2U2l I2U3[ 



5.3. Historical significance of the two-electron problem 

In spite of being the simplest many-body system, the significance of the problem of two 
interacting electrons confined in an external potential cannot be overstated. Historically 
it played a central role in the development of the quantum theory of matter through 
the failure of the Bohr-type semiclassical models to account for the natural He atom. 
Most recently it has influenced the development of several fields like nonlinear physics, 
atomic physics, semiconductror quantum dots, and quantum computing. 

It is instructive to make here a historical detour. Indeed, the failure of Bohr-type 
semiclassical models, based on the orbiting of spatially correlated (antipodal) electrons 
in conjunction with the Bohr-Sommerfeld quantization rule, to yield a reasonable 
estimate of the ground state of the He atom signaled a looming crisis in physics in the 
1920's, which Bohr himself, as well as others, had been keenly aware of, as summarized 
succintly by Sommerfeld: "All attempts made hitherto to solve the problem of the 
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neutral helium atom have proved to be unsuccessful" |205] : see also the 10th chapter 
entitled "It was the Spring of hope, it was the Winter of despair" in the book by Pais 
[206] . the review by Van Vleck ^07\ and the book by Born [208] . 

While, since, numerical solutions of the two-electron Schrodinger equation provided 
a quantitative resolution to the problem, the first successful semiclassical treatment of 
the three-body Coulomb system awaited till 1980 [21)91 [2T0]. 

Furthermore, based on rather general group-theoretical arguments arising from the 
observation of hierarchies with lower symmetry in the excited spectra, and motivated by 
ideas originating in nuclear-physics spectroscopy, it has been discovered in the late 1970's 
and early 1980's that electron correlations in doubly excited He lead to quantization of 
the spectrum much like in a linear triatomic molecule, e-He^^-e. This molecular picture, 
with near rigidity and separability, results in "infinite sequences of vibrational levels, on 
each of which is built an infinite sequence of rotational levels" |202[ 12111 1212j . 



The two previous section [5^1] and section 5^ describing the formation and properties 
of a 2e Wigner molecule in a single QD may be viewed as the culmination of this 
historical background. Interestingly, as in the aforementioned semiclassical treatments, 
the coUinear configuration plays a special role in the molecule-like model, serving perhaps 
as "partial vindication" of the geometry considered originally by Niels Bohr. 



6. Rotating electron molecules in two-dimensional quantum dots under a 
strong magnetic field: The case of the lowest Landau level (uJc/2ujo oo) 

6.1. REM analytic trial wave functions 

In the last ten years, and in particular since 1999 (when it was explicitly demonstrated 
PU] that Wigner crystallization for small systems is related to symmetry breaking at 
the unrestricted Hartree-Fock mean-field level), the number of publications addressing 
the formation and properties of Wigner (or electron) molecules in 2D QDs and 
quantum dot molecules has grown steadily 1201 EH ESI EH EH EQl [Ml HDSl [1281 
[1321 [Ml lini [na [1881 [1891 [1951 ^ A 
consensus has been reached that rotating electron molecules are formed both in zero 
[2Tll23llMlim[5m iT0Tl[T67l[m[TMnMIT95l[2m 
high [2S1 EZl [ini IMl inSl IMl IMl [Ml ESDI E2I] magnetic fields. 

At -B = 0, in spite of considerable differences explored in this report (see next 
paragraph), formation of REMs in quantum dots is driven by the same physical factors 
as Wigner crystallization in infinite 2D media, i.e., when the strength of the interelectron 
repulsion relative to the zero-point kinetic energy (Rw) exceeds a certain critical value, 
electrons spontaneously crystallize around sites forming geometric molecular structures. 
At high magnetic fields, the formation of Wigner molecules may be thought of as 
involving a two-step crystallization process: (I) the localization of electrons results 
from the shrinkage of the orbitals due to the increasing strength of the magnetic field; 
(II) then, even a weak interelectron Coulomb repulsion is able to arrange the localized 
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electrons according to geometric molecular structures (thus this process is independent 
of the value of Rw)- It has been found [211 ISSl 123 I128j that the molecular structures 
at high B coincide with the equilibrium configurations at i? = of classical point 
charges [lEl ESI HH 1222] • 

Due to the small finite number, N, of electrons, however, there are two crucial 
differences between the REM and the bulk Wigner crystal. Namely, (I) the crystalline 
structure is that of the equilibrium 2D configuration of classical point charges, and 
thus consists of nested polygonal rings, in contrast to the well known hexagonal bulk 
crystal; 

(II) in analogy with the case of 3D natural molecules, the Wigner molecules rotate 
as a whole (collective rotations); they behave, however, as highly floppy (non-rigid) 
rotors. 

A most striking observation concerning the REMs is that their formation and 
properties have been established with the help of traditional ab initio many-body 
methods, i.e., exact diagonahzation, [501 lEHl [1311 IMl [HIl EISl |22T] quantum Monte 
Carlo [ZHl HHHl [1951 [Ml |2l6], and the systematic controlled hierarchy [20ll22ll23l 
[211 123 [23 11021 1217] of approximations involving the UHF and subsequent post- 
Hartree-Fock methods. This contrasts with the case of the Jastrow/Laughlin [55] and 
composite-fermion [HSl [FT] wave functions, which were constructed through "intuition- 
based guesswork." 

In spite of its appearance in the middle nineties and its firm foundation in many- 
body theory, however, the REM picture had not, until recently, successfully competed 
with the CF/JL picture; indeed many research papers [HH [2231 [222 1223 12261 [22Z1 [228] 
and books |133j describe the physics of quantum dots in high magnetic fields following 
exclusively notions based on CF/JL functions, as expounded in 1983 (see Ref. [55]) and 
developed in detail in 1995 in Ref. |229] and Ref. |174] . One of the main obstacles for 
more frequent use of the REM picture had been the lack of analytic correlated wave 
functions associated with this picture. This situation, however, changed with the recent 
explicit derivation of such REM wave functions [21]. 

The approach used in Ref. [21] for constructing the analytic REM functions in 
high B consists of two-steps: First the breaking of the rotational symmetry at the 
level of the single-determinantal unrestricted Hartree-Fock approximation yields states 
representing electron molecules. Subsequently the rotation of the electron molecule is 
described through restoration of the circular symmetry via post Hartree-Fock methods, 
and in particular Projection Techniques [18]. The restoration of symmetry goes beyond 
the single determinantal mean-field description and yields multi-determinantal wave 
functions. 

In the zero and low-field cases, the broken symmetry UHF orbitals need to be 
determined numerically, and, in addition, the restoration of the total-spin symmetry 
needs to be considered for unpolarized and partially polarized cases. The formalism 
and mathematical details of this procedure at B = have been elaborated in previous 
sections. 



Symmetry breaking and quantum correlations 



52 



In the case of high magnetic fields, the spins of the electrons are fully polarized. 
Furthermore, one can specifically consider the limit when the confining potential can be 
neglected compared to the confinement induced by the magnetic field, so that the Hilbert 
space is restricted to the lowest Landau level. Then, assuming a symmetric gauge, the 
UHF orbitals can be represented |2H l23Uj by displaced Gaussian analytic functions, 
centered at different positions Zj = Xj + iYj according to the equilibrium configuration 
of N classical point charges |112[ 1113^ 11141 1222j arranged at the vertices of nested regular 
polygons (each Gaussian representing a localized electron). Such displaced Gaussians 
in the lowest Landau level are written as 



u{z,Zj) = {l/y^)exp[-\z - Zj\^/2]exp[-i{xYj -yXj)], (6.1) 

where the phase factor is due to the gauge invariance. z = x + iy, and all lengths 
are in dimensionless units of IbV^ with the magnetic length being Ib = Jhc/eB. 



Note that expression (6.1) is a special case of the more general expression (2.24) for 
a displaced Gaussian which corresponds to situations with smaller magnetic fields when 
the restriction to the lowest Landau level breaks down. The notation z = x + iy is 
associated with positive angular momenta for the single-particle states in the lowest 
Landau level. Ref. used z = x — iy and negative single-particle angular momenta 
in the lowest Landau level. The final expressions for the trial wave functions do not 
depend on these choices. 

Ref. [23] used these analytic orbitals to first construct the broken symmetry UHF 
determinant, and then proceeded to derive analytic expressions for the many- 



Table 1. Comparison of yrast-band energies obtained from REM and EXD 
calculations for iV = 6 electrons in the lowest Landau level, that is in the limit 
_B — > oo. In this limit the external confinement can be neglected and only the 
interaction energy contributes to the yrast-band energies. Energies in units of e^/(K/s). 
For the REM results, the (1,5) polygonal-ring arrangement was considered. The values 



of the fractional filling may be obtained for each L as 


= N{N 


- l)/(2i) 


L 


REM 


EXD 


Error (%) 


L 


REM 


EXD Error (%) 


70 


2.3019 


2.2824 


0.85 


140 


1.6059 


1.6006 


0.33 


75 


2.2207 


2.2018 


0.85 


145 


1.5773 


1.5724 


0.31 


80 


2.1455 


2.1304 


0.71 


150 


1.5502 


1.5455 


0.30 


85 


2.0785 


2.0651 


0.65 


155 


1.5244 


1.5200 


0.29 


90 


2.0174 


2.0054 


0.60 


160 


1.4999 


1.4957 


0.28 


95 


1.9614 


1.9506 


0.55 


165 


1.4765 


1.4726 


0.27 


100 


1.9098 


1.9001 


0.51 


170 


1.4542 


1.4505 


0.26 


105 


1.8622 


1.8533 


0.48 


175 


1.4329 


1.4293 


0.25 


110 


1.8179 


1.8098 


0.45 


180 


1.4125 


1.4091 


0.24 


115 


1.7767 


1.7692 


0.42 


185 


1.3929 


1.3897 


0.23 


120 


1.7382 


1.7312 


0.40 


190 


1.3741 


1.3710 


0.23 


125 


1.7020 


1.6956 


0.38 


195 


1.3561 


1.3531 


0.22 


130 


1.6681 


1.6621 


0.36 


200 


1.3388 


1.3359 


0.21 


135 


1.6361 


1.6305 


0.34 
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body REM wave functions by applying onto an appropriate projection operator 



Vl (see section 2.2.1) that restores the circular symmetry and generates correlated 



wave functions with good total angular momentum L. These REM wave functions 
can be easily written down [21] in second-quantized form for any classical polygonal 



ring arrangement (ni,n2, 



rir 



by following certain simple rules for determining the 



5 Zn J 



where the l/s 



coefficients of the determinants D{li,l2, ■ ■ ■ ,In) = det[zi- , Z2 , 
denote the angular momenta of the individual electrons. 

The REM functions associated with the (0, A^) and (1,A^ — 1) ring arrangements, 
respectively [here (0, N) denotes a regular polygon with vertices, such as an equilateral 
triangle or a regular hexagon, and (1, — 1) is a regular polygon with — 1 vertices 
and one occupied site in its center], are given by 



Zl,Z2, 



X 




1-1, '■2, 



with 



and 



L = Lo + Nm, m = 0,1,2,3, 



(6.2) 



(6.3) 



Zl,Z2, 



,Zn) 




Y]^ sin 

\2<i<j<N 



TT 



A^-l 



(/j Ij 



N 



X D(0,/2,...,/iv)exp(-^2;,z*/2), 



with 



(6.4) 



(6.5) 



L = Lo + {N-l)m, m = 0, 1,2,3, 

where Lq = A^(A^ — l)/2 is the minimum allowed total angular momentum for A^ (fully 
spin polarized) electrons in high magnetic fields. 

Notice that the REM wave functions [Equation (6.2) and Equation (6.4)] vanish 
identically for values of the total angular momenta outside the specific values given by 
the sequences (6.3) and (6.5), respectively; these sequences are termed as magic angular 
momentum sequences. 

We remark that, while the original REM analytic wave function was derived in the 
context of a high magnetic field (that is in the fractional quantum Hall effect regime), 
it is valid for any circumstance where the spectrum consists of a degenerate manifold of 
LLL-like states (even with no magnetic field present). Indeed a wave function having 
the form of the REM wave function discussed by us above has been employed recently 
for graphene quantum dots with a zig-zag boundary condition and in the absence of a 
magnetic field |231j . 

In the remaining of this section, we continue discussing the properties of analytic 
REM wave functions associated with fully spin polarized electrons. However, we mention 
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Figure 17. Total interaction energy from exact-diagonalization calculations as a 
function of the total angular momentum (19 < L < 140) for = 6 electrons in the 
lowest Landau level. The upwards pointing arrows indicate the magic angular momenta 
corresponding to the classically most stable (1,5) polygonal ring arrangement of the 
Wigner molecule. The short downwards pointing arrows indicate successful predictions 
of the composite-fermion model. The medium-size downwards pointing arrows indicate 
predictions of the composite-fermion model that fail to materialize as magic angular 
momenta. The long downward arrows indicate EXD magic angular momenta not 
predicted by the composite-fermion model. Energies in units of / kIb, where k is the 
dielectric constant. 



here that, following the methodology of Ref. [23] for fully spin polarized REMs. Dai et 
al |232] and Shi et al |233] have most recently presented analytic trial wave functions 
for rotating electron molecules with partial spin polarizations. 



6.2. Yrast rotational hand in the lowest Landau level 

As an accuracy test, we compare in table [T] REM and exact-diagonalization results for 
the interaction energies of the yrast band associated with the magic angular momenta 
Lm [see (2.26)] of = 6 electrons in the lowest Landau level. An yraslj|]state is defined 

tt The word yrast is the superlative of the Swedish yr, which means dizzy JJJ. The term yrast is widely 
used in nuclear spectroscopy. 
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as the lowest-energy state for a given angular momentum L. As a result, the yrast 
band represents excitations with purely rotational motion; no other excitations, like 
center-of-mass or vibrational modes, are present. 

As seen from table [T| the REM wave functions offer an excellent approximation to 
the EXD ones, since the relative error of the REM energies is smaller than 0.3%, and 
it decreases steadily for larger L values. Of course, a small difference in the energies 
between approximate and exact-diagonalization results is only one of several tests for 
deciding whether a given trial wave function is a good approximation. As will be 
discussed below, comparison of conditional probability distrubutions is an equally (if 
not more) important test. 

6.3. Inconsistencies of the composite-fermion view for semiconductor quantum dots. 

Before the development of the REM approach, electrons in the lowest Landau level in 
two-dimensional quantum dots were thought as being well approximated by composite 
fermion trial wave functions. However, results obtained with the REM and exact- 
diagonalization calculations led researchers to examine inconsistencies and discrepancies 
of the CP approach in the context of quantum dots. This section focuses on these issues. 



For = 6, figure 17 displays (in four frames) the total interaction energy from 



exact-diagonalization as a function of the total angular momentum L in the range 



19 < L < 140. [The total kinetic energy in the Hamiltonian (4.1), being a constant, 
can be disregarded.] One can immediately observe the appearance of downward cusps, 
implying states of enhanced stability, at certain magic angular momenta. 
For the CF theory, the magic angular momenta can be determined by 

L = L* + mN{N -1) = L* + 2mLo. (6.6) 

Namely, for A^ = 6, if one knows the non-interacting L*'s, the CF magic L's in any 
filling-factor interval l/(2m — 1) > ly > l/(2m -|- 1) [corresponding to the angular- 
momentum interval 15(2m — 1) < L < 15(2m -|- 1)], m = 1,2,3,4,..., can be 
found by adding 2mLo = 30m units of angular momentum to each of the L*'s. To 
obtain the non- interacting L*'s, one needs first to construct f2E \ 1229] the compact 
Slater determinants. Let Nn denote the number of electrons in the nth Landau 
level with X]n=o = N; t is the index of the highest occupied Landau level and 
all the lower Landau levels with n < t are assumed to be occupied. The compact 
determinants are defined as those in which the A^„ electrons occupy contiguously the 
single-particle orbitals (of each nth Landau level) having the lowest angular momenta 
/ = — n, — n -|- 1, . . . , — n -|- A^„ — 1. The compact Slater determinants are usually denoted 
as [A^o, A^i, . . . , Nt\; see Refs. [251 for details. 

The compact determinants and the corresponding non-interacting L*'s for n = 6 
are listed in Table |2l 

There are nine different values of L*'s, and thus the CF theory for A^ = 6 predicts 
that there are always nine magic numbers in any interval 15(2m — 1)<L< 15(2m -|- 1) 
between two consecutive angular momenta of Jastrow/Laughlin states, 15(2m — 1) and 
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15(2m + 1), m = 1, 2, 3, . . . (henceforth we will denote this interval as X^). For example, 



using Table [2j and (6.6), the CF magic numbers in the interval 15 < L < 45 (m = 1) 
are found to be the following nine 

15, 21, 25, 27, 30, 33, 35, 39, 45. (6.7) 

On the other hand, in the interval 105 < L < 135 (m = 4), the CF theory predicts the 
following set of nine magic numbers, 

105, 111, 115, 117, 120, 123, 125, 129, 135. (6.8) 



An inspection of the total-energy-vs-L plots in figure 17 reveals that the 
CF prediction misses the actual magic angular momenta specified by the exact- 
diagonalization calculations as those associated with the downward cusps. It is apparent 
that the number of downward cusps in any interval Im is always different from 9. Indeed, 
there are 10 cusps in Xi [including that at L = 15, not shown in figure Il7| a)], 10 in I2 



[see figure [l7|(b)], 7 in X3 [see figure 17 ^c)], and 7 in X4 [see figure |l7[d)]. In detail, the 
CF theory fails in the following two aspects: (I) There are exact magic numbers that 
are consistently missing from the CF prediction in every interval; with the exception of 
the lowest L = 20, these exact magic numbers (marked by a long downward arrow in 
the figures) are given by L = 10(3m — 1) and L = 10(3m + 1), m = 1, 2, 3, 4, . . .; (II) 
There are CF magic numbers that do not correspond to downward cusps in the EXD 
calculations (marked by medium-size downward arrows in the figures). This happens 
because cusps associated with L's whose difference from Lq is divisible by 6 (but not 
simultaneously by 5) progressively weaken and completely disappear in the intervals 
Xm with m > 3; only cusps with the difference L — Lq divisible by 5 survive. On the 
other hand, the CF model predicts the appearance of four magic numbers with L — Lq 
divisible solely by 6 in every interval X^, at L = 30m =p 9 and 30m =p3, m = l,2,3,... 
The overall extent of the inadequacy of the CF model can be appreciated better by 



Table 2. Compact Slater determinants and associated angular momenta L* ior N — & 
electrons according to the CF presciption. Both L* = —3 and L* — ?> are associated 
with two compact states each, the one with lowest energy being the preferred one. 



Compact state 


L* 


[1,1,1,1,1,1] 


-15 


[2,1,1,1,1] 


-9 


[2,2,1,1] 


-5 


[3,1,1,1] 


-3 


[2,2,2] 


-3 


[3,2,1] 





[4,1,1] 


3 


[3,3] 


3 


[4,2] 


5 


[5,1] 


9 


[6] 


15 
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the fact that there are six false predictions (long and medium-size downward arrows) 
in every interval with m > 3, compared to only five correct ones [small downward 
arrows, see figure 17 'c) and figure [T7|d)]. 

In contrast to the CF model, the magic angular momenta in the REM theory are 
associated with the polygonal ring configurations of N classical point charges. This 
is due to the fact that the enhanced stability of the downward cusps results from the 
coherent collective rotation of the regular-polygon REM structures. Due to symmetry 
requirements, such collective rotation can take place only at magic-angular-momenta 
values. The in-between angular momenta require the excitation of additional degrees 
of freedom (like the center of mass and/or vibrational modes), which raises the total 
energy with respect to the values associated with the magic angular momenta. 

For = 6, the ring configuration of lowest energy is the (1,5), while there exists 
a (0,6) isomer |114l I222j with higher energy. As a result, our exact-diagonalization 
calculations [26] (as well as earlier ones |131t 11321 I171j for lower angular momenta 
L < 70) have found that there exist two sequences of magic angular momenta, a primary 
one (Sp) with L = 15 + 5m [see Equation (6.5)], associated with the most stable (1,5) 
classical molecular configuration, and a secondary one (5*^) with L = 15 + 6m [see 
Equation (6.3)], associated with the metastable (0,6) ring arrangement. Furthermore, 
our calculations (see also Refs. |132[ I171j ) show that the secondary sequence Ss 
contributes only in a narrow range of the lowest angular momenta; in the region of 
higher angular momenta, the primary sequence Sp is the only one that survives and the 
magic numbers exhibit a period of five units of angular momentum. It is interesting to 
note that the initial competition between the primary and secondary sequences, and the 
subsequent prevalence of the primary one, has been seen in other sizes as well |171] i.e., 
= 5, 7, 8. Furthermore, this competition is reflected in the field-induced molecular 
phase transitions associated with broken symmetry UHF solutions in a parabolic QD. 
Indeed, Ref. [53] demonstrated recently that, as a function of increasing B, the UHF 
solutions for = 6 first depict the transformation of the maximum density droplet |126] 
(see definition in section 2.2.1) into the (0,6) molecular configuration; then (at higher 
B) the (1,5) configuration replaces the (0,6) structure as the one having the lower HF 
energy. 

The extensive comparisons in this section lead to the conclusion that the composite- 
fermion model does not explain the systematic trends exhibited by the magic angular 
momenta in 2D quantum dots in high magnetic fields. These trends, however, were 
shown to be a natural consequence of the formation of REMs and their metastable 
isomers. 

These results motivated a reexamination of the original composite-fermion approach 
(the mean-field CF) and led to a reassessment of the significance of the residual 
interaction, neglected in the mean-field CF theory. Initially, it has been reported that 
some CF functions away from the main fractions [e.g., for = 19 and L = 1845 and 
A^ = 19 and L = 3555] may reproduce the aforementioned crystalline patterns 



Subsequently, Jain and coworkers have found that inclusion of the residual 
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interaction is absolutely necessary to account for the full range of inconsistencies of 
the mean- field CF theory |176] . However, this latter development was achieved with 
the trade off of abandoning the fundamental nature of the composite fermion as an 
elementary, independent and weakly interacting quasi-particle. Indeed, the revised |176] 
CF picture amounts to an exact diagonalization method which uses a correlated basis 
set (made out of CF wave functions). 

Another attempt to update the CF theory in order to account for crystallization 
consists of combining the REM analytic wave function ^^^^{zi, Z2, . . . , zn) (see 
section 6.1) with Jastrow prefactors [235], namely one uses a variational wave function 



of the form 

^2p,CFC^-Q^^^_^^.)2p^REM^ (6.9) 

i<j 

with L = L*+pN{N—l) and p serves as an additional variational parameter. Obviously, 
the crystalline patterns in such an approach originate from the REM wave function and 
the Jastrow prefactors simply increase the variational freedom, leading to a numerical 
improvement. Although this approach is a straightforward variational improvement of 
the analytic REM method |24j , it is being referred to |233[ I235j as a composite-fermion 
crystal (CFC). 

More direct variational improvements of the analytic REM wave functions can be 
devised in the spirit of the variation-after-projection method. For example, one can use 
angular-momentum conserving variational parameters in front of the sine coefficients in 
the REM expansion [231] • 

6.4- REM versus Laughlin wave functions: Conditional probability distributions and 
multiplicity of zeroes 



Recent extensive numerical calculations [2M E2] have revealed major disagreements 
between the intrinsic structure of the Jastrow/Laughlin trial wave functions [55] for 
the main fractions u = l/(2m -|- 1) and that of the exact-diagonalization and REM 
wave functions. Indeed, it was found that both EXD and REM wave functions exhibit 
crystalline correlations, while the Jastrow/Laughlin ones are liquid- like as originally 
described in Ref. [SS]. 

To illustrate the differences between the intrinsic structure of the REM and EXD 
states in the lowest Landau level versus the familiar Jastrow/Laughlin ones, we display 



in figure 18 the CPDs for cusp states corresponding to a low filling factor v = 1/5 



and for two different sizes, i.e., for = 6 electrons (L = 75, left column) and N = 7 



electrons {L = 105, right column). In figure 18 the top row depicts the REM case; the 



EXD case is given by the middle row, while the CF case [which reduces to the JL wave 
functions for fractions l/(2p + 1)] are given by the bottom row. 

There are three principal conclusions that can be drawn from an inspection of 



figure 18 (and the many other cases studied in Ref. [2S])- 
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N=6, L=75 (v= 1/5) N=7, L=105 (v= 1/5) 




Figure 18. (Color online) Conditional probability distributions at high B for iV = 6 
electrons and L = 75 (i^ = 1/5, left column) and for iV = 7 electrons and L = 105 
(again v = 1/5, right column). Top row: REM case. Middle row: The case of exact 
diagonalization. Bottom row: The Jastrow/Laughlin case. The exact diagonalization 
and REM wave functions have a pronouned crystalline character, corresponding to 
the (1,5) polygonal configuration of the REM for = 6, and to the (1,6) polygonal 
configuration for N = 7. In contrast, the Jastrow/Laughlin wave functions exhibit a 
characteristic liquid profile that depends smoothly on the number N of electrons. The 
observation point (identified by a solid dot) is located at tq — bASllg for iV = 6 and 
L = 75 and tq = 5.883Zb for TV = 7 and L = 105. The EXD Coulomb interaction 
energies (lowest Landau level) are 2.2018 and 2.9144 e'^/nlg for N = 6,L = 75 
and N — 7, L — 105, respectively. The errors relative to the corresponding exact- 
diagonalization energies and the overlaps of the trial functions with the EXD ones are: 
(I) For iV = 6, L = 75, REM: 0.85%, 0.817; JL: 0.32%, 0.837. (II) For TV = 7, L = 105, 
REM: 0.59%, 0.842; JL: 0.55%, 0.754. 



(I) The character of the exact-diagonahzation states is unmistakably crystaUine 
with the EXD CPDs exhibiting a well developed molecular polygonal configuration 
[(1,5) for = 6 and (1,6) for N = 7, with one electron at the center], in agreement 
with the explicitly crystalline REM case. 

(II) For all the examined instances covering the low fractional fillings 1/9, 1/7, and 
1/5, the Jastrow/Laughlin wave functions fail to capture the intrinsic crystallinity of the 
exact-diagonalization states. In contrast, they represent "liquid" states in agreement 
with an analysis that goes back to the original papers [551 1236] by Laughlin. In 
particular, Ref. |236j investigated the character of the JL states through the use of 
a pair correlation function [usually denoted by g{R)] that determines the probability 
of finding another electron at the absolute relative distance R = |r — ro| from the 
observation point tq. Our anisotropic CPD of Equation (1.1) is of course more general 
(and more difficult to calculate) than the g{R) function of Ref. [236]. However, both 
our P(r, fq) (for = 6 and N = 7 electrons) and the g{R) (for = 1000 electrons, 
and for 1/ = 1/3 and u = 1/5) in Ref. |236] reveal a similar characteristic liquid-like and 
short-range-order behavior for the JL states, eloquently described in Ref. |236j (see p. 
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N=7, L=63 (v=1/3) 




Figure 19. (Color online) CPDs at high B for N = 7 and L = 63 {v = 1/3). 
Top: REM case; Middle: EXD case; Bottom: JL case. Unlike the JL CPD (which is 
liquid) , the CPDs for the exact-diagonalization and REM wave functions exhibit a well 
developed crystalline character [corresponding to the (1,6) polygonal configuration of 
the REM for N = 7 electrons]. The observation point (identified by a solid dot) is 
located at tq = 4.568I_b- 



249 and p. 251). Indeed, we remark that only the first- neighbor electrons on the outer 
rings can be distinguished as separate localized electrons in our CPD plots of the JL 



functions [see figure 18 



(III) For a finite number of electrons, pronounced crystallinity of the exact- 
diagonalization states occurs already at the u = 1/5 value (see figure 18). This finding 
is particularly interesting in light of expectations [234| 1237] (based on comparisons 
between the JL states and the static bulk Wigner crystal) that a liquid- 



to-crystal phase transition may take place only at lower fillings with u < 1/7. 

Of interest also is the case of u = 1/3. Indeed, for this fractional filling, the liquid 
JL function is expected to provide the best approximation, due to very high overlaps 
(better than 0.99) with the exact wave function [58l |59| llGlj . In figure 19, we display 
the CPDs for iV = 7 and L = 63 (z/ = 1/3), and for the three cases of REM, EXD, and 
JL wave functions. Again, even in this most favorable case, the CPD of the JL function 
disagrees with the EXD one, which exhibits clearly a (1,6) crystalline configuration in 
agreement with the REM CPD. 

Similar crystalline correlations at higher fractions were also found for quantum 
dots of larger sizes, e.g., = 8, and = 9 electrons. As illustrative examples for 
these additional sizes [see also the EXD CPD for A^ = 12 electrons in figure 24 below 
(in section 7.3)], we displayed in figure 5 of Ref. [26] the CPDs for A^ = 8 and L = 91 



(1/5 <u = 4/13 < 1/3) and for A^ = 9 and L = 101 ( 1/3 < z/ = 36/101 < 1). Again, 
the CPDs (both for the REM and the EXD wave functions) exhibit a well developed 
crystalline character in accordance with the (1,7) and (2,7) polygonal configurations of 
the REM, appropriate for A^ = 8 and A^ = 9 electrons, respectively. 
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Another area of disagreement between REM and Laughlin wave functions concerns 
the properties of the zero points. In this respect, we recall that the Jastrow/Laughlin 
trial functions for electrons have the form 

$JL(zi,Z2,...,^iv)= I n (^.-^,)''"+Mexp(-f:^,</2). (6.10) 

\i<i<j<N / i=l 

Due to the Jastrow factors (zj — Zj)^'"+^, it is apparent that the Laughlin expressions 



(6.10) (as a function of a given Zi) have N — 1 zero points, each of order 2m + 1, which 
are bound to the positions of the remaining — 1 electrons. In contrast, as discussed in 
Ref. [21], the analytic REM wave functions do not have zeroes with order higher than 
unity. In particular, only A^ — 1 of the REM zeroes are bound to the positions of the 
remaining electrons, while the rest of them are free. Recently, it has been shown through 
extensive numerical studies |239j that the properties of REM zeroes are in agreement 
with the behavior of the zeroes in exact-diagonalization wave functions; this is another 
indication of the superiority of the REM picture compared to the Laughlin theory. 

Before exiting this discussion, we remark about discrepancies of the Laughlin 
quasihole theory in the context of quantum dots. In particular, we recall that the 
Laughlin quasihole, with A^ additional units of angular momentum, has been conjectured 
to be the first excited state. However, LLL exact-diagonalization calculations for A^ 
electrons in a quantum dot have revealed that this is not the case. Instead, the first 
excited state corresponds to an increment in the total angular momentum which varies 
as the number of electrons localized on one of the rings of the rotating electron molecule, 



usually the outermost one; see figure 26 in section 7.4 below. 



7. Rotating electron molecules in two-dimensional quantum dots under a 
strong, but finite external magnetic field (uJc/2ujq > 1) 

7.1. Ground-state energies in medium and high magnetic field 



The general form (2.24) for the displaced Gaussian orbitals [in conjunction with the 



projected REM wave function (2.25)] enables us to calculate REM ground-state energies 
for moderately-high B, when corrections arising from higher Landau levels must be 
taken into consideration. Unlike the lowest-Landau-level case, where the azimuthal 



integration can be carried out analytically, the energies (2.28) (and corresponding CPDs) 



associated with the general REM wave function (2.25 ) require numerical integration over 
the azimuthal angles jq. 

Before proceeding with the presentation of results for A^ > 10, we demonstrate 



the accuracy of the two-step method embodied in Equation (2.25) through comparisons 
with existing exact-diagonalization results for smaller sizes. In figure |20} our REM 
calculations for the ground-state energies as a function of B are compared to EXD 
calculations [171] for A^ = 4 electrons in an external parabolic confinement. The thick 
dotted line (red) represents the broken-symmetry UHF approximation (first step of our 
method), which naturally is a smooth curve lying above the EXD one [solid line (green)]. 
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Figure 20. (Color online) Two-step-method versus exact-diagonalization calculations: 
Ground-state energies for = 4 electrons (referenced to ATioj) as a function of the 
magnetic field B. Thick dashed line (red): broken-symmetry UHF (static electron 
molecule). Solid line (green): EXD (from Ref. |17in . Thick dashed-dotted line (blue): 
REM. Thin dashed line (violet, marked LLL): the commonly used approximate energies 
El^^{B) (see text for details). Thin dotted hue (black): El^^{B) (see text). For 
B < 8 T, the Elf^^{B) and El^l^{B) curves coincide; we have checked that these 
curves approach each other also at larger values of B, outside the plotted range. 
Numbers near the bottom curves denote the value of magic angular momenta [Lm, 
see (2.26)] of the ground state. Corresponding fractional filling factors are specified 
by ^ = N{N — l)/(2L„j). Parameters used: confinement Tlloo = 3.60 meV, dielectric 
constant k — 13.1, effective mass m* — 0.067mp. 



The results obtained after restoration of symmetry [dashed-dotted hne (blue); marked 
as REM] agree very well with the EXD one in the whole range 2T< B < 15T. We recall 
here that, for the parameters of the quantum dot, the electrons form in the intrinsic 
frame of reference a square about the origin of the dot, i.e., a (0,4) configuration, with 
the zero indicating that no electron is located at the center. According to (2.26), Lq = 6, 
and the magic angular momenta are given by = 6-|-4fc, A; = 0, 1, 2, . . . Note that the 
REM energies are slightly lower than the EXD ones in several subranges. According to 
the Rayleigh-Ritz variational theorem, this indicates that the hyperspherical-harmonics 
calculation (equivalent to an exact-diagonalization approach) of |171] did not converge 
fully in these subranges. 



To further evaluate the accuracy of the two-step method, we also display in figure 20 
[thin dashed line (violet)] ground-state energies E^^{B) calculated with the commonly 
used approximate LLL Hamiltonian [HSl 1233 IMI] 

N 2 

HLLL = iV/icD + ;i(cu-^)L+ £ — , (7.1) 

^ j>i=l '^^ij 



where u = \JuJq + ujI/A. The LLL Hamiltonian TYlll reduces to the previously 
introduced Hamiltonian TYlll [see Equation (4.1)] in the limit B ^ oo. Both 
Hamiltonians restrict the many-body wave functions within the lowest Landau level, and 
they both accept the same set of eigenstates as solutions. Indeed the term h{Cj — uJc/2)L 
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is proportional to the total angular momentum, and thus its presence influences only 
the eigenvalues, but not the composition of the eigenstates. Hlll corresponds to a 
situation where the external harmonic confinement is added to TYlll as a perturbation 
(see section II. B in Ref. [53]). As a result, (i) the degeneracy of the single-particle 
levels in the lowest Landau level is lifted and (ii) there is an eigenstate with minimum 
energy (the ground state) at each value of B (expressed through the cyclotron frequency 
Uc). Natur ally, the LLL levels used in the exact diagonalization of Hlll are given by 



expression (4.2), but with A = / = yfi/{m*oo). 

We find that the energies Ef^^{B) tend to substantially overestimate the REM 
(and EXD) energies for lower values of B (e.g., by as much as 5.5% at i? ~ 4 T). On 
the other hand, for higher values of B (> 12 T), the energies E^™{B) tend to agree 
rather well with the REM ones. We stress that the results labelled simply as EXD 
correspond to exact diagonalizations without any restrictions on the Hilbert space, i.e., 
the full Darwin-Fock single-particle spectrum is considered at a given B. 

A behavior similar to E^^^{B) is exhibited also by the E^™{B) ground-state 
energies [which are calcu 
wave functions in section 



6.1 



ated using the Hamiltonian (7.1) and the LLL analytic REM 



with lengths in units of Jfi/{m*uj) instead of IbV^', dotted 



line (black)]. A similar agreement between REM and EXD results, and a similar 



inaccurate behavior of the LLL approximate Hamiltonian ( 7. 1 ) was found by us also 
for = 3 electrons in the range 2T<i?<16T shown in figure 2 of Ref. [53] (the 
exact-diagonalization calculation in this figure was taken from Ref. |166j ). 

In all cases, the total energy of the REM is lower than that of the UHF Slater 



determinant (see, e.g., figure 20). Indeed, a theorem discussed in section 3 of Ref. |241] . 
pertaining to the energies of projected wave functions, guarantees that this lowering of 
energy applies for all values of A^ and B. 

7.2. The case of N = 11 electrons. 



Figure 21 presents the case for the ground-state energies of a quantum dot with A^ = 11 
electrons, which have a nontrivial double-ring configuration (^1,^2). The most stable 
[114] classical configuration is (3,8), for which we have carried UHF (static electron 
molecule) and REM (projected) calculations in the magnetic field range 5 T < i? < 25 T. 
Figure 21 also displays the LLL ground-state energies Effl^{B) [dotted curve (black)]. 



which, as in previous cases, overestimate the ground-state energies for smaller B. The 
approximation Effl^{B), however, can be used to calculate ground-state energies for 
higher values of B. In keeping with the findings for smaller sizes [51] [with (0,A^) or 
(1,A^ — 1) configurations], we found that both the UHF and the REM ground-state 
energies approach, as B ^ 00, the classical equilibrium energy of the (3,8) polygonal 
configuration [i.e., 19.94 meV; 4.865Eo in the units of Ref. [Hi], ^0 = {m*uj^e^ /2KY^''^]■ 



ln analogy with smaller sizes (see, e.g., figure 20 and Ref. [S3]), the REM ground 



state energies in figure 21 exhibit oscillations as a function of B (see in particular 



the inset). These oscillations are associated with magic angular momenta, specified 



Symmetry breaking and quantum correlations 



64 



21.2 




5 10 15 20 25 



B(T) 

Figure 21. (Color online) Ground-state energies for = 11 electrons (per particle, 
referenced to Tioj) as a function of the magnetic field B. Dashed line (red): UHF (static 
electron molecule). Solid line (blue): REM. Dotted line (black): Approximate energies 
E^^l^{B) (see text). Parameters used: confinement UuJn — 3.60 meV, dielectric 
constant k = 13.1, effective mass to* — 0.067TOe. The inset shows a magnification 
of the REM curve in the range 5 T < B < 12 T. 



by the number of electrons on each ring. For = 11 they are given by (2.26), i.e., 
Lm = 55 + 3ki + 8/c2, with the fc^'s being nonnegative integers. As was the case with 
= 9 electrons [S3], an analysis of the actual values taken by the set of indices {/ci, /C2} 
reveals several additional trends that further limit the allowed values of ground-state 
Lm's. In particular, starting with the values {0,0} at i? = 5 T {Lq = 55), the indices 
{^1,^2} reach the values {3,24} at 5 = 25 T (Lm = 256). As seen from table |3| the 
outer index k2 changes faster than the inner index ki. This behavior minimizes the 
total kinetic energy of the independently rotating rings; indeed, the kinetic energy of 
the inner ring (as a function of ki) rises faster than that of the outer ring (as a function 



Table 3. Ground-state magic angular momenta and their decomposition {ki,k2} for 
N = 11 in the nagnetic- field range 5 T < _B < 25 T. The results correspond to the 
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REM (see lower curve in figure 21 1. The parameters used are as in figure 21 
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Figure 22. (Color online) Conditional probability distributions for the REM ground 
state of iV = 11 electrons at i? = 10 T (L = 106). The electrons are arranged 
in a (3,8) structure. The observation point (solid dot) is placed on (left) the outer 
ring at tq = 1.480i?o, and (right) on the inner ring at rg = 0.557i?o. Parameters 
used: confinement huo = 3.60 meV, dielectric constant k = 13.1, effective mass 
m* = 0.067me. Lengths in units of Rq — (2e^/m*Kw§)^/'^. CPDs (vertical axes) 
in arbitrary units. 



of ^2) due to smaller moment of inertia (smaller radius) of the inner ring [see (7.2) . 

In addition to the overestimation of the ground-state energy values for smaller 
magnetic fields (see figure 21 and our discussion above), there are additional 
shortcomings of the lowest-Landau-level approximation pertaining to the ground-state 
ring configurations. In particular, for = 11, we find that according to the LLL 
approximation the ground-state angular momentum immediately after the maximum 
density droplet (Lq = 55) is Lm = 66, i.e., the one associated with the (0, N) vortex-in- 
the-center configuration. This result, erroneously stated in Refs. |242t 1243] as the ground 
state, disagrees with the correct result that includes the full effect of the confinement and 
is listed in table [3} where the ground-state angular momentum immediately following 
the maximum density droplet is = 63. This angular- momentum value corresponds 
to the classicaly most stable (3,8) ring configuration, that is, a configuration with no 
vortex at all (see also the case of = 9 electrons in Ref. [53]). 

Figure |22] displays the REM conditional probability distributions for the ground 
state of = 11 electrons at 5 = 10 T {Lm = 106). The (3,8) ring configuration is 
clearly visible. We note that when the observation point is placed on the outer ring 
(left panel), the CPD reveals the crystalline structure of this ring only; the inner ring 
appears to have a uniform density. To reveal the crystalline structure of the inner ring, 
the observation point must be placed on this ring; then the outer ring appears to be 
uniform in density. This behavior suggests that the two rings rotate independently of 
each other, a property that is explored in the next section to derive an approximate 
quasiclassical expression for the yrast rotational spectra associated with an arbitrary 
number of electrons. 
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7.3. Approximate analytic expression for the yrast-band spectra 



In figure [23| we display ttie CPD for tlie REM wave function of = 17 electrons. This 
case has a nontrivial three-ring structure (1,6,10) |114] which is sufficiently complex to 
allow generalizations for larger numbers of particles. The remarkable floppy character 
(leading to a non-classical, non-rigid rotational inertia, see section VI of Ref. |53j) of 



the REM is illustrated in the CPDs of figure 23 Indeed, as the two CPDs [reflecting 
the choice of taking the observation point [fq in (1.1)] on the outer (left frame) or the 



inner ring (right frame)] reveal, the polygonal electron rings rotate independently of 
each other. Thus, e.g., to an observer located on the inner ring, the outer ring will 
appear as having a uniform density, and vice versa. The wave functions obtained from 
exact diagonalization exhibit also the property of independently rotating rings [see, 
e.g., the = 12 and L = 132 {u = 1/2) case in figure 24 , which is a testimony to the 
ability of the REM wave function to capture the essential physics of a finite number of 
electrons in high B. In particular, the conditional probability distribution displayed in 



figure [24] for exact-diagonalization wave functions exhibits the characteristics expected 
from the CPD evaluated using REM wave functions for the (3,9) configuration and with 



an angular-momentum decomposition into shell contributions [see Equation (2.25) and 
Equation (2.27)] Li = 3 + 3ki and L2 = 63 -|- 9^2 {Li + L2 = Lm', for Lm = 132 the 
angular- momentum decomposition is Li = 6 and L2 = 126). 

In addition to the conditional probabilities, the fioppy-rotor character of the REM 
is revealed in its excited rotational spectrum for a given B. From our microscopic 



calculations based on the wave function in (2.25), we have derived (see below) an 



approximate (denoted as "app"), but analytic and parameter-free, expression [see (7.7) 
below] which reflects directly the nonrigid character of the REM for arbitrary size. 
This expression allows calculation of the energies of REMs for arbitrary N, given the 
corresponding equilibrium configuration of confined classical point charges. 




Figure 23. (Color online) Ground-state conditional probability distributions obtained 
from REM wave functions for the ground state oi N = 17 electrons at B = 10 T 
(L = 228). The electrons are arranged in a (1,6,10) structure. The observation point 
(solid dot) is placed on the outer ring at ro = 1.858i?o (left frame), and on the inner 
ring at tq — 0.969i?o (right frame). The rest of the parameters are: confinement 
fiujQ — 3.6 meV, dielectric constant k = 13.1, effective mass m* = 0.067me. Lengths 
in units of Rq = {2e'^ / {nm* ujq)^/'^ . CPDs (vertical axes) in arbitrary units. 
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Figure 24. (Color online) CPDs for N — 12 electrons and with angular momentum 
L = 132 {v = 1/2) calculated using exact diagonalization in the lowest Landau level. 
The electrons are arranged in a (3,9) structure. The observation point (solid dot) is 
placed on the outer ring at rg — 5.221b (left frame), and on the inner ring at vq — 1.871b 
(right frame). Lengths in units of Ib. CPDs (vertical axes) in arbitrary units. 



We focus on the description of the yrast band at a given B. Motivated by the 
aforementioned nonrigid character of the rotating electron molecule, we consider the 
following kinetic-energy term corresponding to a (ni, . . . , Ug, . . . , n,.) configuration (with 

E^^iN) = E h'Ll/{2JM) - hu;,L/2, (7.2) 

q=l 

where Lg is the partial angular momentum associated with the gth ring about the center 
of the dot and the total angular momentum is L = J2q=i Lq. Jq{aq)) = ngm*a'^ is the 
rotational moment of inertia of each individual ring, i.e., the moment of inertia of 
classical point charges on the qth polygonal ring of radius a^. To obtain the total 
energy, Ef^^^, we include also the term E^^^^N) = J2q=i Jqiflq)^'^ 1'^ due to the effective 
harmonic confinement a) (see Appendix A.l), as well as the interaction energy -E^p, 

^a^pp(^) = E "^^^^ + E E a.). (7.3) 

The first term is the intra-ring Coulomb-repulsion energy of Ug point-like electrons on 
a given ring, with a structure factor 

Sg = E(sm[(j - l)7^/ng])-\ (7.4) 
i=2 

The second term is the inter-ring Coulomb-repulsion energy between rings of uniform 
charge distribution corresponding to the specified numbers of electrons on the polygonal 
rings. The expression fo Vc is 

Vc{ag, a,) = ngn,e^[K{al + a^'/']"' 2Fi[3/4, 1/4; 1; 4a^a2(aJ + al)-% (7.5) 
where 2-^1 is the hypergeometric function. 
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For large L (and/or B), the radii of the rings of the rotating molecule can be found 
by neglecting the interaction term in the total approximate energy, thus minimizing 
only E^lliN) + E^,^^{N). One finds 

aq = XJLq/ng, (7.6) 



with A = / = \Jh/rn*UL)] i.e., the ring radii depend on the partial angular momentum Lg, 

^app' 



reflecting the lack of radial rigidity. Substitution into the above expressions for E^"^ 



E'^pp, and i?^p yields for the total approximate energy the final expression: 

EfSiW = HCo - u^j2)L + j:% + Y.t Vci\i^, (7.7) 

q=l Lq q=l s>q V "9 V "s 

where the constants 

Cv,q = 0.25nJ/25geV(/€A). (7.8) 



For simpler (0, A^) and (1,A^ — 1) ring configurations. Equation (7.7) reduces to the 
expressions reported earlier [51l 1128] . 

The fioppy-rotor character of the REM under strong magnetic field is reflected in 
the absence in (7.7) of a kinetic-energy term proportional to L^. This contrasts with 



the rigid-rotor behavior of an electron molecule at zero magnetic field (see section 5.2 
and Ref. 



7.4- Possible implications for the thermodynamic limit 

While our focus in this section is on the behavior of trial and exact wave functions 
in (finite) quantum dots in high magnetic fields, it is natural to inquire about 
possible implications of our findings to fractional-quantum-Hall-effect systems in the 
thermodynamic limit. 

We recall that appropriate trial wave functions for clean FQHE systems possess 
a good angular momentum L > Lq, a property shared by both the CF/JL and REM 
functions [211 ESI EZl 1236] . We also recall the previous finding [53 I236j that for large 
fractional fillings u > 1/7, the liquid-like (and circularly uniform) Jastrow/Laughlin 
function is in the thermodynamic limit energetically favored compared to the broken- 
symmetry static Wigner crystal (which has no good angular momentum); for u < 1/7, 
the static Wigner crystal becomes lower in energy. This finding was enabled by the 
simple form of the JL functions, which facilitated computations of total energies as a 
function of size for sufficiently large (e.g., N = 1000). 

A main finding of the recent literature on quantum dots is that the exact-numerical- 
diagonalization wave functions of small systems (N < 12) are crystalline in character 
for both low and high fractional fillings. This finding contradicts earlier suggestions 
[551 12291 [236] that, for high z/'s, small systems are accurately described by the liquid-like 
JL wave functions and their descendants, e.g., the composite-fermion ones. Of course, 
for the same high z/'s, our small-size results cannot exclude the possibility that the 
CPDs of the exact solution may exhibit with increasing A^ a transition from crystalline 



Symmetry breaking and quantum correlations 



69 



^ -0.15 

CD 

_E 

>; -0.25 

c 
0) 

c -0.35 
g 



to 
CO 

-0.55 

10 15 20 25 30 

B (T) 

Figure 25. Stabilization energies A_B||-'" = E^-^^ — i?uHF for TV = 6 (dashed curve) 
and N — 7 (solid curve) fully polarized electrons in a parabolic QD as a function of 
B. The troughs associated with the major fractional fillings (1/3, 1/5, and 1/7) and 
the corresponding ground-state angular momenta [L = N{N — l)/(2iy)] are indicated 
with arrows. We have extended the calculations up to _B = 120 T (not shown), and 
verified that AE"!!"" remains negative while its absolute value vanishes as i? ^ oo. 
The choice of parameters is: Tiujq — 3 meV (parabolic confinment), m* = 0.067me 
(electron effective mass), and k — 12.9 (dielectric constant). 




to liquid character, in agreement with the JL function. However, as of now the existence 
of such a transition remains an open theoretical subject. 

For the low fractions, the rotating-electron- molecule theory raises still another 
line of inquiry. Due to the specific form of the REM wave functions, computational 
limitations (in the so-called disk geometry that is natural to quantum dots) prevent 
us at present from making extrapolations of total energies at a given u as N ^ oo. 
Nevertheless, from the general theory of projection operators, one can conclude that 
the REM energies exhibit a different trend compared to the JL ones, whose energies 
were found 1236] to be higher than the static Wigner crystal. Indeed the rotating- 
electron-molecule wave functions remain lower in energy than the corresponding static 
crystalline state for all values of N and z/, even in the thermodynamic limit. This is due 
to an "energy gain" theorem (see Section 3 in Ref. |241] ) stating that at least one of the 
projected states (i.e., the ground state) has an energy lower than that of the original 
broken-symmetry trial function (e.g., the UHF determinant), and this theorem applies 
for any number of electrons and for all values of the magnetic field B. Naturally, 
the REM wave functions will be physically relevant compared to those of the broken- 
symmetry crystal at the thermodynamic limit if the energy gain does not vanish when 
N — > oo] otherwise, one needs to consider the posssibility that the static crystal is the 
relevant physical picture. 

The discussion in the above paragraph may be recapitulated by the following 
question: which state is the relevant one in the thermodynamic limit (A^ oo) - the 
broken-symmetry one (i.e., the static crystal) or the symmetry restored (i.e., rotating 
crystal) state? This question, in the context of bulk broken-symmetry systems, has 
been addressed in the early work of Anderson who concluded that the broken- 
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symmetry state (here the UHF static crystalhne solution) can be safely taken as the 
effective ground state. In arriving at this conclusion Anderson invoked the concept of 
(generalized) rigidity. As a concrete example, one would expect a crystal to behave like 
a macroscopic body, whose Hamiltonian is that of a heavy rigid rotor with a low-energy 
excitation spectrum L'^/2J', the moment of inertia J7 being of order N (macroscopically 
large when — > oo). The low-energy excitation spectrum of this heavy rigid rotor 
above the ground-state (L = 0) is essentially gapless (i.e., continuous). Thus although 
the formal ground state posseses continuous rotational symmetry (i.e., L = 0), "there is 
a manifold of other states, degenerate in the N oo limit, which can be recombined to 
give a very stable wave packet with essentially the nature" [TD] of the broken-symmetry 
state (i.e., the static Wigner crystal in our case). As a consequence of the "macroscopic 
heaviness" as A^ oo, one has: (I) The energy gain due to symmetry restoration (i.e., 
the stabilization energy Ai?|^™ = -E|em ~ -^uhf, see figure 25) vanishes as A^ — > oo, 
and (II) The relaxation of the system from the wave packet state (i.e., the static Wigner 
crystal) to the symmetrized one (i.e., the rotating crystal) becomes exceedingly long. 
This picture underlies Anderson's aforementioned conclusion that in the thermodynamic 
limit the broken-symmetry state may be used as the effective ground state. 

Consequently, in the rest of this section we will focus on issues pertaining to 
the "rigidity" of the rotating electron molecule in high magnetic fields. In particular, 
using our projection method and exact diagonalization, we have demonstrated explicitly 
[50t |5T] that the rigid-rotor picture applies to an A^-electron QD only when B = 0. In 
contrast, in the presence of a high magnetic field, we found [5T 1I521I53] that the electrons 
in the quantum dot do not exhibit global rigidity and therefore cannot be modeled as 
a macroscopic rotating crystal. Instead, a more appropriate model is that of a highly 
non-rigid rotor whose moment of inertia depends strongly on the value of the angular 
momentum L. This behavior originates from the dominance of the magnetic field over 
the Coulomb repulsion. 

The non-rigid rotor at high B has several unique properties: (I) The ground state 
has angular momentum Lgs > 0; (II) While the rotating electron molecule does not 
exhibit global rigidity, it possesses azimuthal rigidity (i.e., all electrons on a given 
ring rotate coherently), with the rings, however, rotating independently of each other. 
Furthermore, the radii of the rings vary for different values of L, unlike the case of a rigid 
rotor; (III) The excitation spectra do not vary as L^; instead they consist of terms that 
vary as aL + Yll=i bq/^fLq (with Ylq=i Lq = L; for the precise values of the constants 
a and b see section 7.3 and Refs. [51] 153]): (IV) The angular momentum values are 



given by the magic values [see section 6.1 L = Lq + Y7q=i kqUq, where (ni, n2, . . . , n^) is 
the polygonal ring arrangement of the static Wigner molecule (with riq the number of 
electrons on the gth ring) and ki < k2 < ■ ■ ■ < kq are nonnegative integers. These magic 
L's are associated with the cusp states which exhibit a relative energy gain with respect 
to neighboring excitations. Thus the low-energy excitation spectrum of the non-rigid 
rotor is not dense and exhibits gaps due to the occurrence of the magic (cusp) states (see 
figure 26). Furthermore, these gaps are reflected in the oscillatory behavior of A£'|g™ 
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Figure 26. Low-energy part of the spectrum of the parabohc QD whose parameters 
are the same as those in figure [25] calculated as a function of the angular momentum 
L through exact diagonalization for N — 7 electrons at a magnetic field B = 18.8 
T. We show here the spectrum in the interval 95 < i < 115 (in the neighborhood 
of = 1/5). The magic angular momentum values corresponding to cusp states are 
marked (99, 105, and 111), and they are seen to be separated from the rest of the 
spectrum. For the given value of B, the global energy minimum (ground state) occurs 
for Lgs = 105, and the gap A to the first excited state {L = 99) is indicated. The 
lowest energies for the different L's (the yrast band) in the plotted range are connected 
by a dashed line, as a guide to the eye. The zero of energy corresponds to Ihib, where 
(D = (wq + w^/4)^/^ and = eB/{m*c). The horizontal arrow denotes the energy 
of the Laughlin quasihole at L=112. It is seen that the Laughlin quasihole is not the 
lowest excited state, as presumed in Ref. [55] . 



(see, e.g., figure 25) as a function of B (or u). 

As increases, more polygonal rings are successively added, and since the 
polygonal rings rotate independently of each other (see, e.g., the case of = 12 in 
figure 24), we expect that the non-rigid-rotor picture remains valid even as A^ oo. 
As a result, it is plausible to conjecture the following properties at high B in the 
thermodynamic limit: (I) the oscillatory character of AE^^^"^ will maintain, yielding 
nonvanishing stabilization energies at the fractional fillings u, and (II) the low-energy 
excitation spectra of the system will still exhibit gaps in the neighborhood of the 
magic angular momenta (see figure 26). Of course, these conjectures need to be 
further supported through numerical calculations for large A^. Nevertheless, the above 
discussion indicates that the question of which state is physically relevant for low 
fractions in the thermodynamic limit at high B - i.e., the broken-symmetry static 
crystal or the symmetrized rotating crystal - remains open, and cannot be answered 
solely following the path of Anderson as described in Ref. [TU]. 

The rotating Wigner crystal has properties characteristic of FQHE states, i.e., it 
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is incompressible (connected to the presence of an excitation gap) and carries a current 
(while the broken-symmetry static crystal is insulating). Thus, we may conjecture 
that a transition at lower fractional fillings from a conducting state with good circular 
symmetry to an insulating Wigner crystal cannot occur spontaneously for clean systems. 
Therefore, it should be possible to observe FQHE-type behavior at low fractional fillings 
in a clean system - a prediction that could explain the observations of Ref. |244j . where 
FQHE behavior has been observed for low fractional fillings typically associated with 
the formation of a static Wigner crystal. In practice, however, impurities and defects 
may influence the properties of the rotating crystal (and its excitations), depending 
on the magnitude of the excitation gap (see, e.g., figure 26). Thus one of the main 
challenges for observation of the fractional quantum Hall effect at such low fillings 
relates to fabrication of high mobility (nearly impurity-free) samples |245j . We remark, 
however, that the stabilization energy and the gap A (see, e.g., figure 26) diminish 
as the magnetic field increases, and as a result the impurities become more efficient 
in influencing the rotating Wigner crystal for the lower fractional fillings (i.e., higher 
angular momenta). 



8. Bosonic molecules in rotating traps: Original results and applications 
8.1. Variational description of rotating boson molecules 

Recent experimental advances in the field of trapped ultracold neutral bosonic gases 
have enabled control of the strength of interatomic interactions over wide ranges 
[851 [86| l87t I246j . from the very weak to the very strong. This control is essential 
for experimental realizations of novel states of matter beyond the well known Bose- 
Einstein condensate [HSl ESI EZl • In this context, the linear ID Tonks- Girardeau regime 
of impenetrable trapped bosons has generated intensive theoretical activity |247[ 1248] 
and several experimental realizations of it have been reported most recently [HSl EH]- 

In this section, we address the properties of strongly-repelling impenetrable bosons 
in rotating ring-shaped or 2D harmonic traps. It has been found that impenetrable 
bosons are "localized" relative to each other [60l 1631 185] and exhibit nontrivial 
intrinsic crystalline correlations j60l [63] . For a small number of bosons, A^, these 
crystalline arrangements are reminiscent of the structures exibited by the well-studied 
rotating electron molecules in quantum dots under high magnetic fields |26l [521 [53] . 
Consequently, we use in the following the term rotating boson molecules. A central result 
of our study is that the point-group symmetries of the intrinsic crystalline structures 
give rise to characteristic regular patterns (see below) in the ground-state spectra and 
associated angular momenta of the RBMs as a function of the rotational frequency for 
neutral bosons (or the magnetic field for charged bosons). 

An unexpected result of our studies is that the rotation of repelling bosons (even 
those interacting weakly) does not necessarily lead to formation of vortices, as is familiar 
from the case of rotating Bose-Einstein condenstates. In particular, for small A^, we will 
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show that the Gross-Pitaevskii energies (including those corresponding to formation of 
vortices) remain always higher compared to the ground-state energies of the RBMs. 
Of course, we expect that the rotating BEG will become the preferred ground state 
for sufficiently large in the case of weakly repelling neutral bosons. We anticipate, 
however, that it will be feasible to test our unexpected results for small by using 
rotating optical lattices, where it is established that a small finite number of atoms can 
be trapped per given site [S7] . 

In a non-rotating trap, it is natural to describe a localized boson (at a position 
Rj) by a simple displaced Gaussian [60]. When the rotation of the trap is considered, 
the Gaussian needs to be modified by a phase factor, determined through the analogy 
between the one-boson Hamiltonian in the rotating frame of reference and the planar 
motion of a charged particle under the infiuence of a perpendicular magnetic field B 
(described in the symmetric gauge). That is, the single-particle wave function of a 
localized boson is 



^j{r) = (/?(r,Rj) = ^^exp 



- ir-(Q X R,) 



with Q = z/A^ and the width of the Gaussian A is a variational parameter; A = 
IbV^ = \/2hc/{eB) = J2h/{muJc) for the case of a perpendicular magnetic field B, 



and A = IqV2 = ^Jh/{mQ,) in the case of a rotating trap with rotational frequency 
(we recall that Uc 2Q, see the Appendix). Note that we consider a 2D trap, so that 
r = {x,y) and R = (X,F). 

The Hamiltonian corresponding to the single-particle kinetic energy is given by 

HKir) = {p-hQx r)V(2m), (8.2) 

for the case of a magnetic field, and by 

HKir) = {p-hQx r)V(2m) - (8.3) 

for the case of a rotating frame of reference jft] 

A toroidal trap with radius tq can be specified by the confining potential 

nr) = ^(r-ro)7^o^ (M 



with Iq = yJi/{muJo) being the characteristic length of the 2D trap. For n ^ 2 and 
^o/fo ~^ this potential approaches the limit of a toroidal trap with zero width, which 
has been considered often in previous theoretical studies (see, e.g., Ref. [249] ) . In the 
following, we consider the case with n = 2, which is more realistic from the experimental 
point of view. In this case, in the limit Tq = 0, one recovers a harmonic trapping 
potential. 



ff The single-particle wave function in (8.1 1 and the many-body projected wave function in (8.51 contain 
contributions from higher Landau levels. These wave functions belong exclusively to the lowest Landau 
level only in the limit when A = V^Ib in the case of a magnetic field, or A = V^hi and fl/ojo = 1 in 
the case of a rotating trap. 
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Figure 27. (Color online) Properties of iV 8 neutral repelling bosons in a rotating 
toroidal trap as a function of the reduced rotational frequency CI/luq. The confining 
potential is given by (8.4| with n — 2 and radius tq = 3Zo, and the interaction-strength 

50. (a) RBM ground-state energies, 

The numbers denote ground-state magic angular 



parameter was chosen as Rs 
shows the range < fi/ojo < 0.3. 
momenta, (b) Energy difference E'^^J — £;UBHF ^ Total angular momenta associated 
with (i) the RBM ground states [thick solid line (showing steps and marked as PRJ); 
online black] and (ii) the UBHF solutions (thin solid line; online red). In the figures, we 
may use the symbol L^, instead of simply L, to denote the 2D total angular momentum. 



To construct an RBM variational many-body wave function describing 
impenetrable bosons in the toroidal trap, we use displaced orbitals v9(r,Rj), 
i = 1,2,...,N [see (8.1)] centered at the vertices of a regular polygon. Then, 
we first construct an unrestricted Bose Hartree-Fock permanent [60l [63] oc 
J2p{im) '^i(jii)¥^2{^i2) ■ ■ ■ V^Af(rijv)- The UBHF permanent breaks the circular symmetry 



of the many-body Hamiltonian. As was discussed in section 3.2, the "symmetry 



dilemma" is resolved through a subsequent "symmetry-restoration" step accomplished 
via projection techniques [231 IMl EHl ED IS21 153] , i-e., we construct a many-body 
wave function with good total angular momentum by applying the projection operator 
Vl = (l/27r) /o^^ deexp[ie{L - L)], so that the final RBM wave function is given by 

\'^n''1) = ^ rde\<l>T''^{e))e'''^. (8.5) 



27r 

l^YP^^i^)) is original UBHF permanent rotated by an azimuthal angle 6. We note 
that, in addition to having good angular momenta, the projected wave function I^I/^r^) 
has also a lower energy than that of [see, e.g. E^^'^ — i^UBHF figure [27|^b). 

The projected ground-state energy is given by 

^PRJ = J^^ h{e)e''''de I J^^ niey- de, (8.6) 
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Figure 28. (Color online) Single-particle densities and CPDs for = 8 bosons in a 
rotating toroidal trap with fl/ujQ = 0.2 and Rs = 50. The remaining trap parameters 
are as in figure 27 (a) Gross-Pitaevskii single-particle density, (b) UBHF single- 
particle density exhibiting breaking of the circular symmetry, (c) RBM single-particle 
density exhibiting circular symmetry, (d) CPD for the RBM wave function [PRJ 



wave function, see Equation (8.5l] revealing the hidden point-group symmetry in the 
intrinsic frame of reference. The observation point is denoted by a white dot. The 
RBM ground-state angular momentum is iz = 16. Lengths in units of ^o- The vertical 
scale is the same for (b), (c), and (d), but difi^erent for (a). 



where h{e) = (^^^^"^(^ = 0)|7^|$^^i^^(^)) and n{d) = = 0)|$]^^"^(^)); the 

latter term ensures proper normahzation. 

The many-body Hamiltonian in the rotating trap is given by 

N N 

n = + V{r,)] + J2v{r,,vj), (8.7) 

i=l i<j 

with the interparticle interaction being given by a contact potential vs{ri, rj) = g6{ri — 
Tj) for neutral bosons and a Coulomb potential vc{ri,rj) = Z^e^/|rj — rj\ for charged 
bosons. The parameter that controls the strength of the interparticle repulsion relative 
to the zero-point kinetic energy is given by Rs = gm/{2'nh^) [201 E3] for a contact 
potential and Rw = Z'^e^ / {TiujqIq) |20l ED] for a Coulomb repulsion. 

For a given value of the dimensionless rotational frequency, VL/ujq, the projection 
yields wave functions and energies for a whole rotational band comprising many 
angular momenta. In the following, we focus on the ground-state wave function (and 
corresponding angular momentum and energy) associated with the lowest energy in the 
band. 

Figure [27|(a) displays the ground-state energy -EpRj of = 8 bosons in a toroidal 
trap as a function of the dimensionless rotational frequency VL/uoq^ with ujq being the 
trap frequency. The prominent features in figure [27|^a) are: (i) the energy diminishes 
as Vt/uiQ increases; this is an effect of the centrifugal force, and (ii) the -EpRj curve 
consists of linear segments, each one associated with a given angular momentum L. 
Most remarkable is the regular variation of the values of L with a constant step of 

c)]. These preferred 

angular momenta L = kN with integer k, are reminiscent of the so called "magic 
angular momenta" familiar from studies of electrons under high-magnetic fields in 2D 
semiconductor quantum dots [261 Ell |53] . 



N units (here = 8) [see inset in figure [27|(a) and figure 27 
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Figure 29. (Color online) Properties of TV = 6 neutral bosons in a rotating harmonic 
trap as a function of the reduced rotational frequency fi/wo- The confining potential 
is given by (8.4 1 with n = 2 and tq — 0, and the interaction-strength parameter was 
chosen as Rs = 50. The intrinsic molecular structure is (1,5). (a) RBM ground-state 
energies, E^^"^ . The inset shows a smaller range. The numbers denote ground-state 
angular momenta, (b) Total angular momenta associated with (i) the RBM ground 
states (thick solid line showing steps; online black) and (ii) the UBHF solutions (thin 
solid line; online red). 



The preferred angular momenta reflect the intrinsic molecular structure of the 
localized impenetrable bosons. We note, that the (0,8) polygonal-ring arrangement 
is obvious in the single-particle density associated with the UBHF permanent [see 
figure [28|^b)]; (0,8) denotes no particles in the inner ring and 8 particles in the 
outer one. After restoration of symmetry, however, the single-particle density is 
circularly symmetric [see the PRJ single-particle density in figure 28 c)] and the 
intrinsic crystallinity becomes "hidden" ; it can, however, be revealed via the conditional 
probability distribution [201 [521 |53l [60] [CPD, see figure gd)]. We note the Gross- 
Pitaevskii single-particle density in figure [28|( a), which is clearly different from the PRJ 
density in figure [28|c). 

The internal structure for charged bosons in a toroidal trap (not shown) is similar 
to that of neutral bosons (figure 28), i.e., a (0,8) ring arrangement, portrayed also in 
the stepwise variation (in steps of 8 units) of the total angular momenta. The internal 
structure is also reflected in the variation of the ground-state total energy as a function of 
the magnetic field. In contrast to the case of neutral bosons, however, the ground-state 
energy curve for charged bosons is not composed of linear segments, but of intersecting 
inverted-parabola-type pieces; this is due to the positive contribution of the Lorentz 
force compared to the negative contribution of the centrifugal force in a rotating trap. 

For RBMs in rotating harmonic traps, the polygonal-ring pattern of localized 
bosons becomes more complex than the simple (0, A^) arrangement that appears 
naturally in a toroidal trap. Indeed, in harmonic traps, one anticipates the emergence 
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Figure 30. (Color online) Properties of GP solutions (thin solid line; online red) 
versus those of RBM wave functions (thick solid line; online green) for iV = 6 neutral 
bosons as a function of the reduced rotational frequency CI/ujq. A harmonic trap is 
considered, and the interaction strength equals Rs = 50. (a) Ground-state energies, 
(b) Associated ground-state angular momenta, (c) GP (BEG) single-particle density 
at fl/uo — 0.65 having 7 vortices with a 6-fold symmetry (thus exhibiting breaking 
of the circular symmetry), (d) RBM single-particle density at fi/wo = 0.65 which 
does not break the circular symmetry, (e) GPD of the RBM at U/luq = 0.65 revealing 
the intrinsic (1,5) crystalline pattern. The white dot denotes the observation point 
ro. Note the dramatic difference in spatial extent between the GP and RBM wave 
functions [compare (c) with (d) and (e). Lengths in units of ^o- The vertical scale is 
the same for (d) and (e), but different for (c). 



of concentric ring structures. For N = 6 neutral bosons in a harmonic trap, we observe 
that, as in the case of a toroidal trap, the ground-state energy as a function of the reduced 
rotational frequency, fi/tuo, [figure [29|(a)] is composed of linear segments, but now the 
corresponding magic angular momenta [figure [29|^b)] vary in steps of — 1 = 5 units. 
This indicates a rotating boson molecule consisting of two polygonal rings; denoted as 
a (1,5) structure, with the inner ring having a single boson and the outer ring five. 

In figure [30|^a), we display the rotating-boson-molecule and mean-field Gross- 
Pitaevskii ground-state energies of = 6 strongly repelling (i.e., Rs = 50) neutral 
bosons in a harmonic trap as a function of the reduced angular frequency of the trap. 
The GP curve (thin solid line; online red) remains well above the RBM curve (thick 
solid line; online green) in the whole range < Q/ujq < 1. The RBM ground-state 
angular momenta exhibit again the periodicity in steps of five units [figure [30|^b)]. As 
expected, the GP total angular momenta are quantized [Lz = (no-vortex) or = 6 
(one central vortex)] only for an initial range < Q/ujq < 0.42. For Q/ujq > 0.42, the 
GP total angular momentum takes non-integer values and ceases to be a good quantum 
number, reflecting the broken-symmetry character of the associated mean field, with 
each kink signaling the appearance of a different vortex pattern of p-fold symmetry 
(p = 1,2,3,4,...) |250] : see an example in figure [30|c). The energetic superiority of 
the RBM wave function over the GP solution demonstrated in figure [30|a) was to 
be expected, since we considered the case of strongly repelling bosons. Unexpectedly, 
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Figure 31. (Color online) Properties of GP solutions (thin solid line; online red) 
versus those of RBM wave functions (thick solid line; online green) for iV = 6 bosons 
as a function of the interaction strength Rg. A harmonic trap is considered, and 
the reduced rotational frequency equals Q/ojo = 0.85. (a) Ground-state energies (b) 
Associated ground-state angular momenta. 



however, for a small number of neutral bosons the energetic advantage of the rotating 
boson molecule persists even for weakly repelling bosons, as illustrated in figure [3l](a). 
Indeed, figure [31]^ a) displays the RBM (thick solid line; online green) and GP (thin solid 
line; online red) ground-state energies for = 6 neutral bosons in a trap rotating with 
Q/ujq = 0.85 as a function of the interaction parameter Rs. The surprising result in 
figure [3T|a) is that the GP energy remains above the RBM curve even for Rs 0. Of 
course the RBM wave function is very close to that of a BEG without vortices when 
Rs —>■ (BEGs without vortices are approximately feasible for small A^). However, 
for small A^, our results show that BEGs with vortices (i.e., for Lz > N) are not the 
preferred many-body ground states; instead, formation of RBMs is favored. Note that 
the energy difference E'^^ — E^^'^ increases rapidly with increasing Rs, reflecting the fact 
that the RBM energies saturate (as is to be expected from general arguments), while the 
GP energies (even with vortices fully accounted for) exhibit an unphysical divergence 
as i?5 ^ oo [figure 31 ^a)]; we have checked this trend up to values of Rs = 100 (not 
shown). Of interest again is the different behavior of the RBM and GP ground state 



angular momenta [figure p31[b)] (see also discussion of figure 30 b)). 

To summarize this section: We have studied the ground-state properties of 
a variational many-body wave function for repelling bosons in rotating traps that 
incorporates correlations beyond the Gross-Pitaevskii mean-field approximations. This 
variational wave function describes rotating boson molecules, i.e., localized bosons 
arranged in polygonal-ring-type patterns in their intrinsic frame of reference. For small 
numbers of neutral bosons, and in particular in the case of GP vortex formation, the 
RBM ground-state energies are lower than those associated with the corresponding 
Gross-Pitaevskii BEG solutions. Given the large differences between the properties 
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of the RBM and BEC wave functions (which become more pronounced for larger 
interaction parameter Rs), and the recently demonstrated ability to experimentally 
control Rs [85l EHl IH3 1246] , we anticipate that our results could be tested in experiments 
involving rotating optical lattices. Detection of rotating boson molecules could be based 
on a variety of approaches, such as the measurement of the spatial extent [contrast the 
RBM and BEC spatial extents in figure [30|(c) - figure [30|e)], or the use of Hanbury 
Brown-Twiss-type experiments |251] to directly detect the intrinsic crystalline structure 
of the RBM. 

8.2. Exact diagonalization for bosons in the lowest Landau level 

Rotating ultracold trapped Bose condensed systems are most commonly discussed in the 
context of formation of vortex lattices, which are solutions to the Gross-Pitaevskii mean- 
field equation [i El [2S21 ESS ESSl ESSl ESZ]. Such vortex lattices have indeed been 
found experimentally for systems containing a large number of bosons |258[ 12591 1260] . 
Nevertheless, several theoretical investigations [671 Ell |69l [TQl [71] of rapidly rotating 
trapped bosonic systems suggested formation of strongly correlated exotic states which 
differ drastically from the aforementioned vortex-lattice states. While experimental 
realizations of such strongly correlated states have not been reported yet, there is already 
a significant effort associated with two-dimenional exact-diagonalization studies of a 
small number of particles (A^) in the lowest Landau level; the LLL restriction corresponds 
to the regime of rapid rotation, where the rotational frequency of the trap Q equals 
the frequency of the confining potential. The large majority [671 [691 COl [71] of such 
exact-diagonalization studies have attempted to establish a close connection between 
rapidly rotating bosonic gases and the physics of electrons under fractional- quant um- 
Hall-effect conditions employing the bosonic version of "quantum-liquid" analytic wave 
functions, such as the Laughlin wave functions, composite-fermion, Moore-Read, and 
Pfaffian functions. 

As described in section [6} the "quantum-liquid" picture for a small number of 
trapped electrons in the FQHE regime has been challenged in a series of extensive studies 
[211 12S1 1121 [SB [S21 [S3| of electrons in 2D quantum dots under high magnetic fields. Such 
studies (both exact-diagonalization and variational) revealed that, at least for finite 
systems, the underlying physical picture governing the behavior of strongly-correlated 
electrons is not that of a "quantum liquid." Instead, the appropriate description is 
in terms of a "quantum crystal," with the localized electrons arranged in polygonal 
concentric rings [21 [2S1 [SB [S21 [S21 II271 II2SI HMj- These "crystalline" states lack 
[521 [53] the familiar rigidity of a classical extended crystal, and are better described 
[211 [2ni 1121 [SB [S21 [S3] as rotating electron (or Wigner) molecules. 

Motivated by the discovery in the case of electrons of REMs at high B (and from 
the fact that Wigner molecules form also at zero magnetic field [2Q1 [2S1 [IB [SOl I167[I188] ) 
some theoretical studies have most recently shown that analogous molecular patterns 
of localized bosons do form in the case of a small number of particles inside a static 
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or rotating harmonic trap |l3l IMl ESI I177[ 1178] . In analogy with the electron case, 
the bosonic molecular structures can be referred to [63] as rotating boson molecules; 
a description of RBMs via a variational wave function built from symmetry-breaking 
displaced Gaussian orbitals with subsequent restoration of the rotational symmetry was 



presented in Refs. [131 EDI ES] and reviewed in section 8.1 



In a recent paper, Baksmaty et al [72] used exact diagonalization in the lowest 
Landau level to investigate the formation and properties of RBMs focusing on a larger 
number of particles than previously studied, in particular for sizes where multiple-ring 
formation can be expected based on our knowledge of the case of 2D electrons in high 
B. A finite number of particles (A^ < 11) at both low (z/ < 1/2) and high {u > 1/2) 
filling fractions u = N{N — 1)/2L (where L = C/h is the quantum number associated 
with the total angular momentum C) was studied and both the cases of a long-range 
(Coulomb) and a short-range (5-function) repulsive interaction were investigated. In 
this section, we report some main results from Ref. [72] . 

As in the case of electrons in 2D quantum dots, we probe the crystalline nature of the 
bosonic ground states by calculating the full anisotropic two-point correlation function 
P(r,ro) [see Equation (1.1)] associated with the exact wavefunction \E'(ri, r2, . . . , Ttv). 
The quantity P{r, fq) is proportional to the probability of finding a boson at r given 
that there is another boson at the observation point fq, and it is often referred to as the 
conditional probability distribution (section 1.5). A main finding of our studies is that 
consideration solely of the CPDs is not sufficient for the boson case at high fractional 
fillings u > 1/2; in this case, one needs to calculate even higher-order correlation 
functions, e.g., the full A^-point correlation function defined as the modulus square 
of the full many-body EXD wave function, i.e., 

P(r; n, r2, . . . , rjv-i) = |^(r; ri, ra, . . . , r,v-i)|', (8.8) 

where one fixes the positions of A^ — 1 particles and inquiries about the (conditional) 
probability of finding the A^th particle at any position r. 

The investigations in this section are also motivated by recent experimental 
developments, e.g., the realization of trapped ultracold gas assembhes featuring bosons 
interacting via a long-range dipole-dipole interaction [2611 1262] . We expect the 
results presented in this section to be directly relevant to systems with a two-body 
repulsion intermediate between the Coulomb and the delta potentials. Additionally, we 
note the appearance of promising experimental techniques for measuring higher-order 
correlations in ultra-cold gases employing an atomic Hanbury Brown-Twiss scheme [251] 
or shot-noise interferometry [2631 1264] . Experimental realization of few-boson rotating 
systems can be anticipated in the near future as a result of increasing sophistication 
of experiments involving periodic optical lattices co-rotating with the gas, which are 
capable of holding a few atoms in each site. A natural first step in the study of such 
systems is the analysis of the physical properties of a few particles confined in a rotating 
trap with open boundary conditions (i.e., conservation of the total angular momentum 
L). 
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Figure 32. Ground-state angular momenta, Lg 



for = 6 bosons in a rapidly 
as a function of the 



rotating trap [described by the LLL Hamiltonian in (8.9 1] 
rotational frequency f2 expressed in units of o^o- The bosons interact via a Coulombic 
repulsion (left) and a delta repulsion (right), and the many-body Hilbert space is 
restricted to the lowest Landau level. The angular momentum associated with the first 
bosonic Laughlin state occurs dX L = 30, i.e., at N{N — 1). The value of c = 0.2?ia;oA 
for the Coulomb case (left) and the value oi g — 2nT%ijjQh? / N for the case of a delta 
repulsion (right); the many-body wave functions do not depend on these choices. In 
the delta-interaction case, the values of the angular momenta terminate with the value 
i = 30 (the Laughlin value) at fi/wo = 1- In contrast, in the Coulomb-interaction 
case (left), the values of the ground-state angular momenta do not terminate, but 
diverge as fi/wg 1. Note the stepwise variation of the values of the ground-state 
angular momenta in both cases, indicating the presence of an intrinsic point-group 
symmetry associated with the (0,6) and (1,5) polygonal-ring structure of a rotating 
Boson molecule. 



The main results of Ref. ^72j can be summarized as follows: Similar to the well- 
established (see section |6] and section [?]) emergence of rotating electron molecules in 
quantum dots, rotating boson molecules form in rotating harmonic traps as well. The 
RBMs are also organized in concentric polygonal rings that rotate independently of 
each other, and the polygonal rings correspond to classical equilibrium configurations 
and/or their low-energy isomers. Furthermore, the degree of crystallinity increases 
gradually with larger angular momenta L's (smaller filling fractions z/'s), as was the 
trend [26l |52| 153] for the REMs and as was observed also for v < 1/2 in another 
study |178] for rotating bosons in the lowest Landau level with smaller and single- 
ring structures. We finally note that the crystalline character of the RBMs appears to 
depend only weakly on the range of the repelling interaction, for both the low (see also 
Ref. |178] ) and high (unlike Ref. |177] ) fractional fillings. 

In studies of 2D quantum dots, CPDs were used some time ago in Refs. 
[50l 11281 1131] . For probing the intrinsic molecular structure in the case of ultracold 
bosons in 2D traps, however, they were introduced only recently by Romanovsky et 
al [on]- The importance of using CPDs as a probe can hardly be underestimated. 
Indeed, while exact-diagonalization calculations for bosons in the lowest Landau level 
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Figure 33. (a) Single-particle densities [n{r); left column], (b) CPDs [P(r,ro)] in 
3D plots (middle column), and (c) CPDs in contour plots (right column), portraying 
the strengthening of the crystalline RBM structure for = 6 bosons interacting via 
a repulsive Coulomb interaction as the filling fraction is reduced. The white dots 
in the CPD plots indicate the reference point Tq. We note in particular the gradual 
enhancement of the peak at the center of the plots, and the growth of the radius of 
the outer ring; the latter reflects the nonrigid-rotor nature of the RBMs (in analogy 
with the findings of Ref. [52] regarding the properties of rotating electron molecules). 
The cases of = 1/4 and v — 1/8 exhibit a clear (1,5) crystalline arrangement, while 
the case oi v — 1/2 (first Laughlin state) is intermediate between a (1, 5) and a (0, 6) 
pattern (see text for details). Lengths in units of A. The vertical scales are in arbitrary 
units, which however do not change for the panels within the same column (a), (b), or 
(c). 



have been reported earlier [671 EHl EHl [TOl [7T] , the analysis in these studies did not include 
calculations of the CPDs, and consequently formation of rotating boson molecules was 
not recognized. 

8.2.1. The case of N = 6 bosons in the lowest Landau level. As a specific example of 



the points discussed above in section 8.2, we present here results for = 6 bosons in 



the lowest Landau level. For additional cases (e.g., = 9 and A^ = 11), see Ref. 



In analogy with the magnetic- field Hamiltonian of Equation (7.1), the many-body 
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Hamiltonian for bosons in a rotating trap is reduced in the lowest Landau level to 
the expression 

TV 

n^LL = Nhujo + Hcoo - n)L + Y.v{vi,vj), (8.9) 

i<j 

where uq specifies the 2D harmonic trap and Q denotes the rotational frequency. The 
interparticle interaction is given by a contact potential vs{ri, rj) = g5{vi~Yj) for neutral 
bosons and a Coulomb potential Vciji, r^) = c/\vi — Vj\ for charged bosons. 

Since T^lll is rotationally invariant, i.e., [TYlll^-^] = 0? eigenstates must 
also be eigenstates of the total angular momentum with eigenvalue hL. For a given 
rotational frequency f2, the eigenstate with lowest energy is the ground state; we denote 
the corresponding angular momentum as Lgg. 

We proceed to describe the EXD results for = 6 particles interacting via a 
Coulomb repulsion by referring to figure 32, where we plot Lgs against the angular 
frequency Vt of the rotating trap. A main result from all our calculations is that Lgs 
increases in characteristic (larger than unity) steps that take only a few integer values, 
i.e., for A^ = 6 the variations of Lgs are in steps of 5 or 6. In keeping with previous 
work on electrons [Ml EHl HSJ EH |52l [53] at high and very recently on bosons in 
rotating traps jUJ |60l [631 1177[ I178j . we explain these magic- angular-momenta patterns 
(i.e., for A^ = 6, Lgs = Lq + 5/c or Lgg = Lq + 6fc, with Lq = 0) as manifestation of 
an intrinsic point-group symmetry associated with the many-body wave function. This 
point-group symmetry emerges from the formation of RBMs, i.e., from the localization 
of the bosons at the vertices of concentric regular polygonal rings; it dictates that the 
angular momentum of a purely rotational state can only take values Lgs = Lq + I]i kifii, 
where rii is the number of localized particles on the ith polygonal ring. [We remind the 
reader that for spin-polarized electrons in the lowest Landau level, the corresponding 
value is Lq = N{N — l)/2.] Thus for A^ = 6 bosons, the series Lgs = 5k is associated 
with an (1, 5) polygonal ring structure, while the series Lgs = 6k relates to an (0, 6) 
arrangement of particles. It is interesting to note that in classical calculations |114j for 
A^ = 6 particles in a harmonic 2D trap, the (1, 5) arrangement is found to be the global 
energy minimum, while the (0, 6) structure is the lowest metastable isomer. This fact 
is apparently reflected in the smaller weight of the Lgs = Qk series compared to the 
Lgs = 5k series, and the gradual disappearance of the former with increasing L. 

Magic values dominate also the ground state angular momenta of neutral bosons 
(delta repulsion) in rotating traps, as shown for A^ = 6 bosons in the right panel 



of figure 32 Although the corresponding fi-ranges along the horizontal axis may be 
different compared to the Coulomb case, the appearance of only the two series 5k 
and Qk is remarkable — pointing to the formation of RBMs with similar (1, 5) and 
(0,6) structures in the case of a delta interaction as well (see also Refs. |177[ 1178] ). 
An important difference, however, is that for the delta interaction both series end at 
Q/ujq = 1 with the value L = A^(A^ — 1) = 30 (for A^ = 6 the bosonic Laughlin value 
at z/ = 1/2), while for the Coulomb interaction this L value is reached for Q/ujq < 1 — 
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allowing for an infinite set of magic angular momenta [larger than N{N — 1)] to develop 
as Q/uq ^ 1. 

Beyond the analysis of the ground-state spectra as a function of Q, the intrinsic 
crystalline point-group structure can be revealed by an inspection of the CPDs [and to 
a much lesser extent by an inspection of single-particle densities]. Because the EXD 
many-body wave function is an eigenstate of the total angular momentum, the single- 
particle densities are circularly symmetric and can only reveal the presence of concentric 
rings through oscillations in the radial direction. The localization of bosons within the 
same ring can only be revealed via the azimuthal variations of the anisotropic CPD 
[Equation (1.1)]. One of our findings is that for a given the crystalline features in the 
CPDs develop slowly as L increases (or u decreases). 

For u < 1/2, we find that the crystalline features are well developed for all sizes 
studied by us. In figure 33, we present some concrete examples of CPDs from exact- 
diagonalization calculations associated with the ground-states of = 6 bosons in a 
rotating trap interacting via a repulsive Coulomb potential. In particular, we present 
the CPDs for Lgs = 30 (bosonic Laughlin for u = 1/2), 60, and 120; these angular 



momenta are associated with ground states at specific i7-ranges [see figure 32 . All 



three of these angular momenta are divisible by both 5 and 6. However, only the 
Lgs = 30 CPD (figure 33 top row) has a structure that is intermediate between the 
(1, 5) and the (0, 6) polygonal-ring arrangements. The two other CPDs, associated with 
the higher Lgs = 60 and Lgs = 120 exhibit clearly only the (1,5) structure, illustrating 
our statement above that the quantum-mechanical CPDs conform to the structure of 
the most stable arrangement [i.e., the (1, 5) for A^ = 6] of classical point-like charges as 
the fractional filling decreases. 

However, for u > 1/2, the azimuthal variations may not be visible in the CPDs, 
in spite of the characteristic step-like ground-state spectra [see figure 32 for A^ = 6 
bosons]. This paradox is resolved when one considers higher-order correlations, and 
in particular A^-point correlations [see Equation (8.8)]. In figure 34 and figure 35, we 
plot the A^-point correlation functions for A^ = 6 bosons and Lgg = 15 for both the 
Coulomb interaction and 5-repulsion, respectively, and we compare them against the 
corresponding CPDs. The value of 15 is divisible by 5, and one expects this state to 
be associated with a (1,5) molecular configuration. It is apparent that the CPDs fail 
to portray such fivefold azimuthal pattern. The (1,5) pattern, however, is clear in the 
A^-point correlations (middle and right panels). One has two choices for choosing the 
positions of the first five particles (white dots), i.e., one choice places one white dot at 
the center and the other choice places all five white dots on the vertices of a regular 
pentagon. For both choices, as shown by the contour lines in the figures, the position of 
maximum probability for the sixth boson coincides with the point that completes the 
(1, 5) configuration [see the black dots in the middle and right panels]. 

Note that the differences in the CPDs and A^-point correlation functions between 
the Coulomb and the 5-repulsion are rather minimal. 
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Figure 34. Contour plots of the CPD (a) and A''-point correlation function (b) and 

(c) for iV = 6 bosons with Lgs = 15 interacting via a Coulomb repulsion. The white 
squares indicate the positions of the fixed particles. The black square in (b) and (c) 
indicates the position of the 6th particle according to the classical (1,5) molecular 
configuration. Note the different arrangements of the five fixed particles, i.e.. (b) one 
fixed particle at the center and (c) no fixed particle at the center. Note also that the 
CPD in (a) fails to reveal the (1,5) pattern, which, however, is clearly seen in the 
AT-point correlation functions in both (b) and (c) . Lengths in units of A. The vertical 
scales are arbitrary, but the same in (b) and (c). 



DELTA INTERACTION N=6 1^=15 
(a) CPD y (b) N-POINTl y (c) N-POINT2 




-3 3 -3 3 -3 3 



Figure 35. Contour plots of the CPD (a) and iV-point correlation function (b) and 
(c) for A'' = 6 bosons with Lgg = 15 interacting via a J-repulsion. The white squares 
indicate the positions of the fixed particles. The black square in (b) and (c) indicates 
the position of the 6th particle according to the classical (1,5) molecular configuration. 
Note the different arrangements of the five fixed particles, i.e., (b) one fixed particle at 
the center and (c) no fixed particle at the center. Note also that the CPD in (a) fails 
to reveal the (1,5) pattern, which, however, is clearly seen in the A'^-point correlation 
functions in both (b) and (c). Lengths in units of A. The vertical scales are axbitraxy, 
but the same in (b) and (c). 
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9. Summciry 

This report reviewed the physics of strong correlations in two-dimensional small finite- 
size condensed- matter systems, such as electrons in quantum dots and repelling bosons 
in harmonic traps. It was shown that strong correlations in such systems relate to 
the appearance of symmetry breaking at the mean-field level of description. Particular 
attention was given to the similarities of symmetry breaking in these systems despite 
the different interparticle interactions (Coulombic repulsion in quantum dots versus a 
contact potential for neutral bososns in harmonic traps). 

The universal aspects of symmetry breaking in small systems (including nuclei and 
molecules in quantum chemistry) have been exploited to develop a two-step method 
of symmetry breaking at the unrestricted Hartree-Fock level and subsequent symmetry 
restoration via post Hartree-Fock projection techniques. In conjunction with exact- 
diagonalization calculations, the two-step method was used to describe a vast range of 
strongly-correlated phenomena associated with particle localization and formation of 
crystalline (molecular) structures of electrons in quantum dots and bosons in harmonic 
traps. Due to their finite size, these crystalline structures are different from the familiar 
rigid crystals of extended systems; they rather resemble and exhibit similarities with 
the natural 3D molecules (e.g., ro- vibrational spectra). 

It was shown that strongly-correlated phenomena emerging from symmetry 
breaking include: 

(I) Chemical bonding, dissociation, and entanglement in quantum dot molecules 
and in electron molecular dimers formed within a single anisotropic quantum dot, with 
potential technological applications to solid-state quantum-computing devices. 

(II) Electron crystallization, with localization on the vertices of concentric polygonal 
rings, and formation of rotating electron molecules in circular quantum dots. At zero 
magnetic field, the REMs can approach the limit of a rigid rotor; at high magnetic field, 
the REMs are highly floppy, with the rings rotating independently of each other. 

(III) In the lowest Landau level, the rotating electron molecules are described by 
parameter-free analytic many-body wave functions, which are an alternative to the 
composite- fermion and Jastrow/Laughlin approaches, offering a new point of view of 
the fractional quantum Hall regime in quantum dots (with possible implications for the 
thermodynamic limit). 

(IV) Crystalline phases of strongly repelling bosons. In the case of rotating traps 
and in analogy with the REMs, such repelling bosons form rotating boson molecules, 
which are energetically favored compared to the Gross-Pitaevkii solutions, in particular 
in the regime of vortex formation. 
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Figure Al. (Color online) The Darwin-Fock single-particle energy levels of a 2D 
harmonic oscillator under the effect of a perpendicular magnetic field _B as a function 
of 77 = (jJc/wo, where ujc is the cyclotron frequency and ujo is the frequency specifying 
the 2D harmonic confinement. A specific color (online) indicates orbitals with the 
same number of radial nodes. 



Materials Science. 



Appendix A. 

In this Appendix, we briefly review the single-particle wave functions and associated 
energy spectra of a two-dimensional circular harmonic oscillator under the influence 
of a perpendicular magnetic fleld B (relevant to the case of quantum dots) or under 
rotation with angular frequency Q (relevant to the case of trapped atomic gases in 
rotating harmonic traps). These single-particle wave functions and associated spectra 
are known as the Darwin-Fock states and energy levels, after the names of the authors 
of two original papers [162^ 1163] on this subject. 



Appendix A.l. 
magnetic field 



Two-dimensional isotropic oscillator in a perpendicular 



In this case, the Hamiltonian (for an electron of mass m*) is given by: 
I , e . .0 1 



H 



2m* 



IP 



(A.l) 



where r = (x, y) and ujq is the frequency of the oscillator. In the symmetric gauge, the 
vectror potential is given by A = (B x r)/2, and the Hamiltonian ( A.l ) can be rewritten 
in the form 



H 



ujJ H — ^ ^ r , 

2m* 2 2 



(A.2) 
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where / = —ih{xd/dy — yd/dx) is the angular momentum operator of the electron (in 
the z direction), Uc = eB/{m*c) is the cyclotron frequency, and uj = + u'^/i is the 
effective-confinement frequency. 

The eigenf unctions of the Hamiltonian (A. 2) have the same functional form as 
those of a 2D harmonic oscillator at B=0, but with an effective frequency uj, i.e., in 
polar coordinates 

<Pn,i{P,0) = Ar„,,pl'le-^'/2e^'^Ll|l(p2), (A.3) 



with p = r/l; the characteristic length / is given by / = yh/{m*Co). In (A.3), the index 
n denotes the number of nodes in the radial direction, and I (without any subscript 
or tilde) denotes the angular-momentum quantum numbers; the i^iJ''s are associated 
Laguerre polynomials. 

The single-particle energy spectrum corresponding to the Hamiltonian (A. 2) is 
plotted in figure |A1[ the associated eigenenergies are given by 



nl 



(2n+ |/| + 1)W1 



(A.4) 



with T) = cJc/i^o- 

In the limit of Uc/ (2ujq 
energy spectrum in Equation (A.4) reduces to that of the celebrated Landau levels, i.e., 

(A.5) 



oo, one can neglect the external confinement, and the 



E 



M 



where A4 = n + (|/| — l)/2 is the Landau- level index. 

We remark that the Landau levels are infinitely degenerate. The lowest Landau level 
= contains all nodeless levels [n = 0) with arbitrary positive angular momentum 
/ > 0. 



Appendix A. 2. 



Two-dimensional rotating harmonic oscillator 



In the case of a rotating isotropic oscillator, instead of the expression (A. 2), one has the 
following single-particle Hamiltonian: 

1 

2m 



H 



ni 



2 2 

2 ° ' 



(A.6) 



where the mass of the particle (e.g., a bosonic or fermionic atom) is denoted by m; Q 
denotes the rotational frequency. 



From a comparison of the second terms in (A. 2) and (A.6), one derives the 
correspondence Q uJc/2. 

We note that, unlike the application of a perpendicular magnetic field, the rotation 
does not generate an effective confinement different from the original external one 



[compare the third terms between (A. 2) and (A.6)]. As a result, the eigenfunctions 
of the Hamiltonian (A.6) are given by the expressions (A.3), but with p = t/Iq where 



the characteristic length /g 




Figure A2. (Color online) The Darwin-Fock single-particle energy levels of a 2D 
harmonic oscillator rotating with angular frequency f2 as a function of = U,/ujq, 
where is the frequency specifying the 2D harmonic confinement. A specific color 
(online) indicates orbitals with the same number of radial nodes. 



The single-particle energy spectrum corresponding to the Hamiltonian (A. 6) is 
plotted in figure |A2| and the associated eigenenergies are given by 

(2n+|/| + l)-/r7, (A.7) 



with 7] = VL/uoq. 

For VL/ujq = 1, the energy spectrum in Equation (A.7) reduces to that of the 
corresponding Landau levels, i.e., 

Em = 2hiJo{M + ^), (A.8) 

where A4 = n + (|/| — /)/2 is the Landau-level index. 

As was the case with the perpendicular magnetic field, the Landau levels are in- 
finitely degenerate, and the lowest Landau level = contains all nodeless levels 
{n = 0) with arbitrary positive angular momentum / > 0. However, unlike the magnetic- 
field case where huc depends on B, the energy gap between the Landau levels in the 
case of rotation is independent of Q and equals 2huJo. 
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